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ABSTRACT 


In  this  thesis  the  performance  of  a  step  quantum  well  infrared  photodetector,  designed 
by  Kevin  Lantz  (June  2002)  and  experimentally  studied  by  Michael  Touse  (September  2003) 
and  Yeo  Hwee  Tiong  (December  2004),  was  simulated  in  Matlab  using  the  transfer  matrix 
method.  The  results,  obtained  by  the  Matlab  program,  are  compared  with  the  experimental 
results,  in  an  attempt  to  make  inferences  about  the  optimum  way  of  designing  QWIP  detec¬ 
tors. 

Simulation  of  the  above  implies  numerical  solution  of  the  Schrodinger  equation,  us¬ 
ing  algorithms  and  methods,  which  give  accurate  results.  In  our  approach,  the  transfer  matrix 
method  (TMM)  was  used  with  exponentials  and  Airy  functions  to  represent  the  solutions  to 
Schrodinger  equation  under  zero  and  non-zero  bias,  respectively.  The  calculated  results  were 
compared  with  the  experimental  data  and  found  to  provide  a  good  agreement  which  validated 
the  accuracy  of  the  model  employed. 

In  the  final  section  of  the  thesis  we  examine  and  simulate  in  Matlab  the  application  of 
the  extended  Kalman  filtering  (EKF)  to  an  infrared  photodetector  as  a  target  tracking  mecha¬ 
nism  to  both  maneuvering  and  non-maneuvering  targets.  When  we  used  one  sensor  for  track¬ 
ing,  the  results  were  reliable  provided  that  the  target  did  not  maneuver.  In  the  case  of  a  ma¬ 
neuvering  target  the  results  were  significantly  improved  when  we  used  both  sensors  for 
tracking. 
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I.  INTRODUCTION 


A.  QUANTUM  WELL  INFRARED  PHOTODETECTORS  (QWIPS)-BRIEF 

REVIEW 

In  the  recent  years  there  has  been  a  considerable  increase  in  research  activities 
towards  the  development  of  quantum  well  based  infrared  sensors  [1-4].  The  use  of  quan¬ 
tum  wells  made  of  larger  bandgap  materials,  can  overcome  the  difficulties  associated 
with  smaller  bandgap  materials,  typically  used  in  conventional  infrared  detectors  [1]. 

The  smaller  bandgap  materials  are  relatively  unstable,  and  difficult  to  process,  resulting 
in  non-uniformities  when  focal  plane  arrays  are  made  [1],  The  infrared  detectors  have 
many  commercial  and  military  applications. 


The  operation  of  quantum  well  infrared  photodetectors  (QWIPs)  is  based  on  the 
excitation  of  bound  electrons  in  a  quantum  well  by  infrared,  as  illustrated  in  Figure  1.1a. 
The  quantum  wells  are  formed  by  sandwiching  a  small  bandgap  material  between  large 
bandgap  materials,  as  shown  in  Figure  1.1b.  It  is  possible  to  use  the  quantum  well 
formed  either  in  the  conduction  band  or  valence  band  to  fabricate  QWIPs.  In  the  case  of 
conventional  infrared  detectors,  the  electrons  in  the  valence  band  are  excited  across  the 
bandgap,  while  in  QWIP  detectors  it  is  necessary  to  dope  the  quantum  well  to  populate 
the  ground  state  using  either  n-type  (for  the  conduction  band  well),  or  p-type  (for  the  va¬ 
lence  band  well). 


Figure  1.1  (a)  excitation  of  bound  electrons  due  to  incident  photon  energy  fico  .  (b) 

Quantum  well  formation  by  sandwiching  the  appropriate  bandgap  material  (From  Ref. 

[5]-) 


1 


In  the  design  of  quantum  wells,  the  parameters,  such  as  its  dimension,  and  the 
composition  of  each  structure  material  used,  allow  us  to  manipulate  the  behavior  of  the 
photoexcited  carriers,  so  that  they  escape  from  the  potential  wells  (quantum  leak)  and  are 
collected  as  current  by  the  application  of  an  external  bias,  as  shown  in  Figure  1.2. 


Figure  1.2  Conduction  band  diagram  of  a  QIWP  in  an  external  electric  field.  Excited 
electrons  escape  from  the  ground  state  creating  a  photocurrent  (From  Ref.  [4].) 

The  layers  of  QWIP  detectors  are  primarily  grown  by  molecular  beam  epitaxy 
(MBE),  which  gives  thickness  control  in  atomic  scale  [2,  4].  For  achieving  good  material 
quality,  the  lattice  constant  of  different  materials  used  in  QWIPs  should  be  nearly  the 
same.  This  avoids  the  dislocations  generated  by  the  lattice  mismatch. 

The  usual  wavelength  coverage  of  QWIPs  is  between  3  and  20  microns  [1,8].  By 
controlling  the  barrier  height  and  well  width,  it  is  possible  to  adjust  the  energy  level  sepa¬ 
ration  and,  hence,  the  wavelength  dependence  of  the  detector  response.  For  example,  the 
longer  wavelength  response  can  be  achieved  using  shallow  quantum  wells. 

QWIPs  can  be  categorized  according  to  the  electron  resulting  state  in  four  types: 
bound  to  bound  (B-B),  bound  to  quasibound  (B-QB),  bound  to  continuum  (B-C),  and 
bound  to  miniband  (B-M),  as  schematically  shown  in  Figure  1.3. 
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Figure  1.3  Energy  band  diagram  showing  B-B,  B-QB,  B-C,  and  B-M  QWIP  struc¬ 
tures  (From  Ref.  [4].) 

The  main  factor  which  limits  the  QWIP  performance  is  the  dark  current  (the  cur¬ 
rent  that  flows  through  a  biased  detector  in  the  dark  with  no  photons  impinging  on  it)  , 
which  plagues  all  infrared  detectors.  At  temperatures  above  45  K,  the  dark  current  of  the 
QWIP  is  entirely  dominated  by  classic  thermionic  emission  of  ground  state  electrons  di¬ 
rectly  out  of  the  well  into  the  energy  continuum.  There  are  several  theoretical  techniques 
of  minimizing  it  [3],  most  importantly  the  positioning  of  the  first  excited  state  to  align 
with  the  barrier. 

Recent  studies  indicate  that  multicolor  operation  of  QWIPs  is  possible  by  simply 
stacking  different  types  of  quantum  wells.  A  four-color  QWIP  was  reported  [5],  based  on 
stacked InGaA si  AlGaAs  and GaAs!  AlGaAs  multi-quantum  wells.  In  addition,  a  great 
deal  of  QWIP  focal  plane  arrays  have  been  analyzed  and  demonstrated  by  many  groups 
[6-8],  showing  excellent  performance  in  LWIR  atmospheric  window. 

B.  PURPOSE  OF  THIS  THESIS 

The  purpose  of  this  thesis  was  to  simulate  in  Matlab  the  performance  of  symmet¬ 
ric  and  asymmetric  quantum  well  structures  under  both  zero  bias  and  non-zero  bias.  It 
uses  a  design  made  in  earlier  thesis  by  Kevin  R.  Lantz  [25]  and  experimental  measure¬ 
ments  made  by  Michael  P.  Touse  [26]  and  Yeo  Hwee  Tiong  [27].  The  results  obtained  by 
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the  MATLAB  program  are  compared  with  the  experimental  results,  in  an  attempt  to  make 
inferences  about  the  optimum  way  of  designing  QWIP  detectors.  Simulation  of  the  above 
implies  numerical  solution  of  the  Schrodinger  equation,  using  algorithms  and  methods 
which  give  accurate  results.  Various  methods  have  been  presented  to  solve  the 
Schrodinger  equation  numerically  including  the  WKB  approximation,  the  Monte  Carlo 
method  and  the  finite  element  method,  which  all  basically  deal  with  solving  second  order 
differential  equations.  In  our  approach,  the  transfer  matrix  method  will  be  used  [12,  16] 
with  exponential  and  Airy  functions  to  represent  the  solutions  to  Schrodinger  equation 
under  zero  and  non-zero  bias,  respectively. 

In  the  final  section  of  the  thesis,  we  will  examine  and  simulate  in  Matlab  the  ap¬ 
plication  of  the  extended  Kalman  filtering  (EKF)  to  an  infrared  photodetector  as  a  target 
tracking  mechanism,  to  both  maneuvering  and  non-maneuvering  targets.  Using  a  recur¬ 
sive  type  algorithm  such  as  EKF,  we  will  examine  the  bearing  only  tracking  (BOT)  prob¬ 
lem  classically.  The  main  drawback  of  BOT  is  that  the  research  is  closed,  meaning  that 
the  publication  available  is  limited. 

C.  MILITARY  RELEVANCE 

QWIPs  have  a  great  deal  of  military  and  security  applications.  The  most  common 
are  night  vision,  remote  sensing  and  satellite  IR  imaging,  FLIR,  surveillance,  targeting  in 
a  variety  of  terrains  (air  sea  ground)  and  weapon  delivery  [4,7].  They  are  also  widely 
used  in  search  and  rescue  situations,  mine  detection,  missile  seeker  fonnulation,  smart 
munitions,  weapon  sights,  preventive  maintenance,  non-destructive  testing,  medical  im¬ 
aging  and  surveillance  [4],  Finally,  the  bearing-only-tracking  application  examined  in 
the  last  chapter  of  the  thesis,  is  considered  a  modern  real-world  problem  featuring  the  ad¬ 
vantage  of  passive  target  acquisition  and  surveillance.  Since  modern  military  needs  are 

oriented  in  systems  which  can  track  targets  without  being  noticed,  the  use  of  a  modern 
QWIP  sensor  combined  with  a  classic  tracking  algorithm,  such  as  EKF,  could  result  to  a 
reliable  surveillance  system. 

D.  THESIS  OUTLINE 

The  present  thesis  starts  with  a  general  discussion  of  the  necessity  of  QWIPs,  fol¬ 
lowed  by  a  brief  review  of  how  these  devices  are  formulated  and  where  they  find  applica¬ 
tions  of  military  interest.  Next,  the  reader  is  led,  through  the  mathematical  background,  to 
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the  simulation  program,  and  the  comparison  to  the  laboratory  case.  Finally  a  “real  world” 
application  is  investigated  through  BOT. 

In  particular,  the  outline  of  the  present  thesis  is  as  follows: 

Chapter  II  discusses  the  Schrodinger  equation  and  the  formulation  of  the  transfer 
matrix  method.  The  Schrodinger  equation  is  applied  in  solving  the  unbiased  and  biased 
quantum  well,  and  the  transfer  matrix  method  is  applied  in  order  to  expand  the  solution  in 
more  than  one  region  of  the  structure. 

Chapter  III  provides  information  about  the  constructed  Matlab  code,  and  the  labo¬ 
ratory  semiconductor,  while  simulations  are  presented  to  show  the  QWIP  behavior  to 
several  conditions  of  bias.  Ideas  such  as  probability  current  and  bound-to-continuum  tran¬ 
sitions  are  presented,  and  finally  the  comparison  to  the  laboratory  case  takes  place. 

Chapter  IV  is  focused  into  the  BOT  problem.  The  used  initialization  of  the  state 
vector  is  presented,  as  long  as  the  EKF  algorithm,  with  emphasis  into  bearing  tracking. 
Finally,  we  present  and  discuss  several  cases  of  sensor-target  geometry. 

Chapter  V  is  a  summary  of  the  results,  followed  by  suggestions  for  probable  fu¬ 
ture  research. 

Appendix  A  includes  the  Matlab  coding  used  to  model  the  QWIP. 

Appendix  B  includes  the  Matlab  coding  used  to  model  the  BOT  problem. 

Appendix  C  is  an  overview  of  Airy  functions  and  their  asymptotic  expressions 
used  to  model  the  QWIP  in  the  low  bias  cases. 

Finally,  Appendix  D  includes  the  complete  set  of  resulting  graphs  of  the  BOT 
problem. 
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II.  MATHEMATICAL  AND  PHYSICAL  BACKGROUND 


A.  INTRODUCTION 

The  infrared  photodetectors  examined  in  this  thesis  are  quantum  mechanical  de¬ 
vices.  Modeling  this  type  of  device  requires  solution  of  the  key  quantum  mechanical 
equation,  Schrodinger’ s  equation.  Given  the  potential  and  boundary  conditions,  we  de¬ 
termine  the  shape  and  the  temporal  evolvement  of  the  wave  function  by  solving  this  sec¬ 
ond-order,  partial  differential  equation.  Due  to  the  one  dimensional  nature  of  quantum 
well  structures  used  in  QWIPs,  it  is  sufficient  to  solve  the  Schrodinger  equation  along  the 
growth  direction  of  the  layers. 

Over  the  years,  various  methods  have  been  used  for  numerical  solution  of  the 
Schrodinger  equation,  namely,  the  WKB  approximation  [9],  the  variational  calculation 
method  [10,  1 1],  the  Monte  Carlo  method  [12],  the  finite  element  method  (FEM)  [13], 
and  the  transfer  matrix  method  (TMM)  [14].  In  the  present  work,  we  employed  the  trans¬ 
fer  matrix  method  (TMM)  due  to  its  flexibility  in  analyzing  complex  quantum  well  struc¬ 
tures  under  external  bias. 

B.  SCHRODINGER  EQUATION 

The  one  dimensional  time-dependent  Schrodinger  equation  for  an  electron  can  be 
written  as: 

-- +  U(x)^(x,  t )  =  m  ¥(x,  t),  (2.1) 

2  me  ox  ot 

where  fi  is  the  reduced  Plank’s  constant,  me  is  the  electron  mass,  x  is  the  coordinate 
(usually  taken  along  the  direction  of  the  growth),  'P  is  the  wavefunction  and  U  is  the  po¬ 
tential  energy. 

Since  the  quantum  well  potential  is  independent  of  time,  it  is  possible  to  separate 
the  spatial  and  temporal  dependencies  with  the  spatial  dependence  satisfying  the  time- 
independent  Schrodinger  equation: 

-^I^T(x)  +  U(x)W(x)  =  fPF(x) .  (2.2) 

2  m  ox 


7 


c. 


EFFECTIVE  MASS  MODEL 


The  periodic  potential  of  the  atoms  in  semiconductor  materials  influences  the  pro¬ 
pagation  of  an  electron  inside.  In  order  to  include  the  effect  of  the  periodic  potential,  the 
electron  mass  m  is  usually  replaced  by  an  effective  mass  m  ,  and  treats  the  electron  as  a 
free  particle  [1].  The  effective  mass m  is  related  to  the  band  structure  by: 


m 


fi2 


d  E  /  dk 


2  ’ 


(2.3) 


where  fi  is  the  Plank’s  constant,  E  is  the  conduction  band  energy  and  k  the  wavevector. 


The  QWIPs  are  usually  made  using  potential  wells  in  the  conduction  band  well 
due  to  smaller  effective  mass  which  gives  higher  absorption  [1].  For  the  materials  used 
in  this  thesis  ( InGaAs  ,AlGaAs  )  the  effective  mass  of  electrons  for  various  composi¬ 
tions  is  shown  in  Figure  2. 1 


Figure  2.1:  Effective  mass  as  a  function  ofxfor  Alfia^As  and/ntG2/,  xAs  . 
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D.  TRANSFER  MATRIX  METHOD  (TMM) 

1.  General 

The  transfer  matrix  method  is  a  relatively  simple  approach  for  solving  Schrodin- 
ger  equation  for  quantum  well  structures.  It  is  generally  used  for  a  number  of  problems 
mostly  involving  second-order,  partial  differential  equations  (e.g.,  propagation  character¬ 
istics  in  planar  waveguides,  spontaneous  emission  in  DFB  lasers)  [12-14,  16]. 

Comparing  this  method  with  the  other  methods  usually  used  for  solving  the  same 
problems,  TMM  has  the  advantage  both  of  simplicity  and  usage  in  low  processing  ability 
computers.  Its  simplicity  mainly  involves  multiplication  of  2x2  matrices  which  are 
straightforward  to  implement  on  a  computer.  We  will  examine  two  cases  using  the  TMM, 
quantum  well  structures  with  and  without  the  application  of  an  external  electric  field. 

2.  Quantum  Well  Structure  without  Applied  Electric  Field 

We  assume  an  arbitrary  multilayered  quantum  well  structure  as  shown  in  Figure 

2.2 


Figure  2.2:  Multilayered  quantum  well  structure  with  zero  applied  bias. 


The  time  independent  Schrodinger  equation  for  each  square  quantum  well  region 
of  the  multilayered  quantum  well  structure  reads: 


n2  a2vF,.(x) 

,  .  a  +F.(x)TC(x)  =  EW  (x) .  (2.4) 

2m  dXj  1  1  1 
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where  fi  is  the  reduced  Plank’s  constant,  m *  is  the  electron  effective  mass  in  the  region  j , 

x  is  the  coordinate  (usually  taken  along  the  direction  of  the  growth),  'P  is  the  wave-  func¬ 
tion,  Vj  is  the  potential  energy  in  region  j ,  and  E  stands  for  the  associated  electron  en¬ 
ergy.  The  general  solution  of  Equation  (2.17)  in  each  region  is  a  superposition  of  the  left 
and  the  right  traveling  waves  and  is  given  by: 


¥  ,(*)=  C/  '  ~  D.e  ik 


(2.5) 


where  kj  is  implemented  from  the  separation  of  variables  and  represents  the  wavenumber 


and  C j,Dj  are  arbitrary  constants.  The  value  of  kj  depends  on  the  region  effective  mass 
m*  and  the  potential  height  V,  in  the  region  j  as: 


h 


(2.6) 


At  each  boundary  the  probability  current  should  be  continuous  which  implies: 


_L  d_ 

m  *  dx 


L  J  mj+1  dx  ■  J 


(2.7) 

(2.8) 


where  x .  is  the  coordinate  at  the  boundary  between  j  and  j  + 1  layers. 


Applying  the  solution  of  Equation  (2.4),  given  by  Equation  (2.5),  to  Equations 
(2.7)  and  (2.8),  we  can  find  after  some  algebra: 


k;  ,m:  +  k.m. 


k,  ,m,  -k.m, 


V  -•  ,  ,  T  /V  .rn.  1  :K  x  _:U  r  IV  •  ,  i  III-  IV  -III  •  ,  i  _jhx  r 

Q  _  J  +  X  J  j  J  +  l  Q  j  j  Q  lKj+'xj  Q  |  J  +  l  J  j  ./  +  1  c  lKjxj  e  ixj+ix!  Y) 


2kj+tmJ 


2kj+lmJ 


kMm,  —km,,  ikx  ik  ,  kMmj+kmi+,  _ikx  +ik  x 
U  _  1  +  1  }  ]  }+i  £  J  i  C  1  C  |  7  +  1  7  1  1  +  1  g  j  £  j+'  J  f) 


2kj+lmJ 


J 


2kj+\mj 


J 


(2.9) 


(2.10) 


The  results  in  Equations  (2.9)  and  (2.10)  can  be  expressed  in  a  matrix  form  as 

(2.11) 


=  M . 

~Ci~ 

dj4 

J 

Di_ 
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where  M  j  is  given  by  as 


f  ,  *  \ 

f  i  *  5 

\ 

l+kJmJ+ 1 

h' 

+ 

1 

1  kjmJ+ 1 

e-i(kj+kj+1)Xj 

l  kMm*j) 

l  kMm)\ 

(  J  *  ^ 

1  kjmJ+ 1 

e'(kj  +kj4i)Xj 

l  +  kjmj+ 1 

e  ,a  *  i|r 

k,,m* 

(2.12) 


Repeatedly  applying  Equation  (2. 1 1),  it  is  possible  to  relate  the  coefficients  if  the 
outennost  layers  as: 


r  c 

JV 

rci 

^  mn 

rci 

=  M  .M .  ... 

. MM, 

= 

1Dn] 

J  J  1 

La: 

lm2i 

m22) 

LaJ 

In  the  case  of  bound  states,  the  wavefunction  must  decrease  exponentially  to  zero 
at  the  outennost  boundaries.  This  implies  that  exponential  growing  coefficients  in  Equa¬ 
tion  (2.5),  namely  ( DN,Cl ),  should  equate  to  zero.  Under  this  condition  Equation  (2.13) 
reduces  to: 


~c 

'-'N 

'mn 

mt2' 

~  0 

0 

Mu 

m22 J 

A 

(2.14) 


It  can  be  easily  seen  from  Equation  (2.14)  that,  to  satisfy  the  boundary  conditions, 
we  demand  the  matrix  element  m21  =  0 .  Since  m22  depends  only  on  energy  E ,  the  bound 

states  can  be  found  by  plotting  m22  ( E )  and  identifying  the  zero  crossings. 

3.  Quantum  Well  Structure  under  Applied  Electric  Field 
a.  General 

During  the  QWIP  operation,  it  is  necessary  to  apply  external  bias  to  extract 
photoexcited  electrons.  Under  an  external  electric  field,  the  quantum  well  potential  tilts, 
as  illustrated  in  Figure  2.3.  The  amount  of  the  tilt  depends  on  the  strength  and  the  direc¬ 
tion  of  the  applied  field.  The  bias  alters  the  energy  levels  in  the  quantum  well  structure, 
and  hence,  the  photoresponse.  In  the  section  following,  the  effect  of  an  external  electric 
field  on  QWIP  operation  is  analyzed  by  solving  Schrodinger  equation  with  a  linear  poten¬ 
tial  in  addition  to  the  quantum  well  potential. 
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Figure  2.3  Multilayered  quantum  well  structure  with  forward  applied  bias  (left),  and 

reverse  applied  bias  (right). 


b.  Solution  of  the  Schrodinger  Equation  in  a  Linear  Potential 

We  first  consider  the  solution  of  Schrodinger  equation  in  a  linear  potential 
as  shown  in  Figure  2.4. 


Using  Schrodinger’ s  Equation  (2.2)  and  substituting  the  value  of 
U (x)  =  -ax  we  can  get  to  the  following  form: 


h1  82 

2m  fix' 


-  ¥  (x)  -  axY  (x)  =  E'P  (x) . 


(2.15) 


It  is  possible  to  convert  Equation  (2.15)  to  a  standard  differential  equation 
by  using  the  following  change  of  variables. 

x — y 


x  =  py+y  <=>  y  = 


P 


(2.16) 
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Using  the  parameters  as  defined  in  Equation  (2.20),  the  Equation  (2.19) 
can  be  transformed  into  the  Airy’s  differential  equation: 

-fYW-yW  =  0,  (2.21) 

dy 

which  has  two  linearly  independent  solutions  and  are  given  in  terms  of  Airy  functions  of 
the  first  and  the  second  kind  ( Ai  and  Bi ).  The  general  solution  can  be  written  as  a  linear 
combination  of  the  two  Airy  function  multiplied  by  two  arbitrary  constants  C  and  D  [16]: 

H/(T)=  CAi(y)+DBi(y).  (2.22) 

Figures  (2.5)  and  (2.6)  show  the  Airy  functions  of  the  first  kind,  Ai ,  and 
the  second  kind,  Bi ,  as  a  function  of  y  for  both  positive  and  negative  directions. 
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Figure  2.5:  Airy  functions  of  the  first  kind. 


Figure  2.6:  Airy  functions  of  the  second  kind. 
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In  Equation  (2.16)  by  replacing  the  values  of  y  and  J3  we  get: 


y  = 


X  + 


a 


r  n2  > 

1/3 

\2ma  j 

(2.38) 


where  v  =  0  or  x  =  -E/a  corresponding  to  the  classical  turning  point  as  shown  in  Figure 


Figure  2.7 :  Classical  turning  point  at x  = -E/a  and  v  =  E  . 


Using  the  solutions  of  the  Airy  equation  given  by  (2.22)  and  replacing  the 
variable  v  ,  the  solution  in  terms  of  the  original  coordinate  x,  the  energy  E ,  the  Con¬ 
stantsa  and  fi ,  the  particle  mass  m  and  the  constants  C  and  D  can  be  written  as: 


E 

E 

x  +  — 

x  -\ — 

a 

a 

'F(x)  =  CAi 

(  n2  ] 

1/3 

+DBi 

f  n2  ] 

1/3 

y2ma  y 

v2  ma  y 

(2.39) 


The  constants  C  and  D  are  determined  by  the  boundary  conditions  and  for 
the  above  potential  D  =  0  since  Bi  is  exponentially  growing  in  the  barrier  region  or  nega¬ 
tive  x  direction. 
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c.  Square  Quantum  Well  Under  an  Electric  Field 

Generally,  when  an  electric  field  is  applied  to  a  quantum  well  structure  as 
schematically  illustrated  in  Figure  2.8,  the  profile  of  the  potential  will  be  changed  [1,  12]. 


Figure  2.8:  Single  quantum  well  structure  with  forward  applied  bias. 

This  change  is  given  by  the  equation: 

V(x,F)  =  V(x,0)-exF,  (2.40) 

where  V (x,  0)  is  the  potential  profile  of  the  quantum  well,  F  is  the  applied  electric  field 
in  V/m ,  e  is  the  electron  charge  and  x  is  the  associated  spatial  coordinate. 

By  substitution  of  Equation  (2.40)  into  the  Schrodinger  equation  we  arrive 
at  the  following  fonnula: 

Pi2  <32vF(x)  ,  x 

~  ,  a  2  +  (Vj-exF  (x)T7.(x)  =  i?y.(x).  (2.41) 

2  m.(x)  dXj  v  ’ 

Since  the  potential  energy  VJ  is  a  constant  within  the  jth  layer  of  the  quan¬ 
tum  well,  the  solution  to  Equation  (2.41)  can  be  written  as  a  linear  combination  of  Airy 
functions  [16]  as: 

'F(x)  =  CiAi(pi(x))+  DjBifaixj) ,  (2.42) 

and  p.  (x)  is  defined  as 
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-Fx-rjj 

PF)  =  (  ,  y'3  ’ 

h-F2 

2m* e 

V  i 

where  m*  is  the  effective  mass  associated  with  the  jth  layer  and  ///  is  given  by: 

(E-VA 


where  E  is  the  energy  of  the  electron. 


(2.43) 


(2.44) 


The  boundary  conditions  given  by  Equations  (2.7)  and  (2.8)  are  still  appli¬ 
cable  to  the  Airy  functions  solution,  given  by  Equation  (2.42).  An  arbitrary  coordinate 
system  was  chosen  to  be  at  the  upper  left  end  of  the  potential  well,  for  simplicity  pur¬ 
poses,  as  we  illustrated  in  Figure  2.8. 

The  boundary  conditions  between  two  neighbor  layers  j  and  j  + 1  of  the 
quantum  well,  give  the  following  relations  between  the  coefficients: 


C/+1  -  2..I/  (//;..)/>’/(//.  )j  +  />  |/h( /;,)/>>/(//)- 


(2.45) 


C.+1  |2  ,l/(/f),l/  (//,  ) -  ,l/(//;),l/(//;.;)!  +  D/+]  j2;.l/(//.  )/>7  (/>'.,)-  Ai'(f)  )/>’/(//  ,)! 


(2.46) 


where  /? +1  and  (5j  are  defined  as: 


/?/+i  =  (~Xj  -77/+l  ) , 


2m*  ef 


2m*  ef 


and  F  is  the  ratio  between  the  effective  masses  of  the  two  layers. 


(2.47) 


(2.48) 
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Equations  (2.45)  and  (2.46)  can  be  written  in  matrix  form  as: 


vAy 


n 


Ai(/l:..)/JHfi  )~;:Ai(fi  )  Bi(P:..)Bi(P  )-EBHP,)Bi(P:..) 

ZjAi(Pj)Ai(pj+ 1)  -  Ai\pj)Ai(pj+l)  ^Ai^ppBi  (Pj+l)~  Ai  (P j)Bi{p] /+1) 


fC  ^ 


(2.49) 


KDm  ) 


A  repeated  application  of  boundary  conditions,  similar  to  that  of  the  unbi¬ 
ased  case,  the  coefficients  of  the  outermost  layers  can  be  related  as: 

(2.50) 


...M .  M, 

~cN  1 

Pnu 

mn 

CN 

La: 

J  1  J 

LavJ 

lm2l 

m22y 

[dn\ 

In  the  case  of  bound  states,  the  wavefunctions  in  the  outermost  layers 
should  decay,  which  implies  that  the  coefficients  that  make  them  grow  ( Dl  and  Cv ) 
should  equate  to  zero.  Thus,  Equation  (2.50)  reduces  to: 


1 

_ 1 

Pnu  ml27 

1 

o 

1 _ 

1 

O 

i _ 

\m2\  m22  J 

1 

Q 

_ 1 

(2.51) 


This  implies  that, 


m22  (E)  =  0 . 


(2.52) 


which  corresponds  to  quantized  energy  states  E  inside  the  well  under  the  electric  field. 

Figure  2.9  represents  a  single  AI03Ga07As  quantum  well  of  40  Angstroms 

width  where  the  applied  field  is2x  107  V/m  .  We  represent  the  potential  of  the  well  using 
a  blue  line,  the  wavefunction  with  the  black  line,  and  the  energy  level  with  the  red  line. 
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X  (Angstroms) 

Figure  2.9:  Single  quantum  well  structure  with  large  forward  applied  bias 

(2xl07  V/m). 

For  example,  the  bound  state  energies  in  the  quantum  well  were  obtained 
by  plotting  m12(E )  as  a  function  of  energy,  as  shown  in  Figure  2.10,  and  identifying  the 

energies  where m22(E)  =  0  .  For  the  parameters  used,  there  was  one  bound  state  in  the 

quantum  well,  at  energy  of  about  0.08  9  eV  .  The  extent  of  the  wavefunction  beyond  the 
right  side  of  the  barrier  indicates  tunneling  through  the  barrier  due  its  finite  width.  This 
is  responsible  for  leakage  current  in  quantum  well  detectors  at  low  temperatures. 
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Graph  of  matrix  element  a22  Vs  Energy  E 
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Figure  2.10:  Zeros  of  the  matrix  element  m21  detennines  the  energy  characteristic 

values. 


E.  SUMMARY 

In  this  chapter,  we  developed  the  necessary  tools  to  comprehend  the  solution  of 
the  Schrodinger  equation  in  a  quantum  well  structure,  with  and  without  the  presence  of 
electric  field.  The  use  of  transfer  matrix  method  with  Airy  functions  enabled  the  analysis 
of  quantum  well  structures  under  external  bias.  In  the  next  chapter,  response  of  a  QWIP 
made  from  InxGat  xAs  /  GaAs  step  quantum  well  will  be  simulated  and  compared  with 

available  experimental  data.  Ideas  such  as  continuous  states  and  “normalization”  of  them, 
Fermi’s  golden  rule  for  optical  transitions  and  probability  current  will  be  developed. 
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III.  EMULATION  ANALYSIS  AND  DATA 


A.  INTRODUCTION 

The  simulation  program  developed  in  Chapter  II  was  employed  to  analyze  a 
QWIP  detector  designed  by  Kevin  Lantz  [25]  and  experimentally  studied  by  Michael 
Touse  [26],  and  Yeo  Hwee  Tiong  [27].  The  comparison  between  the  analysis  and  ex¬ 
perimental  results  is  focused  on  the  infrared  absorption  and  detector  responsivity. 

B.  DESCRIPTION  OF  QWIP  SAMPLE 

The  fabricated  step  QWIP  consists  of  a  multiple  quantum  well  structure.  The 
layer  specification  of  the  complete  detector  is  presented  in  Table  1. 


Substance 

Mole  % 

Thickness 

(A) 

Si  Doping 
Concentration 
(cmA-3) 

GaAs 

5000 

1E+18 

GaAs 

300 

InGaAs 

10 

40 

InGaAs 

30 

40 

1E+18 

GaAs 

500 

GaAs 

8000 

1E+18 

GaAs 

6.35E+06 

Table  1.  Specifications  of  the  test  QWIP  (From  [26].) 


To  increase  the  active  layer  thickness  of  the  detector,  25  step  InxGax_xAs  quantum 

wells  were  employed.  A  single  period  of  the  step  quantum  well  is  schematically  shown  in 
Figure  3.1. The  width  of  the  well,  and  the  step,  are  40  Angstroms  each,  while  a  300- 
Angstrom  barrier  of  GaAs  was  placed  between,  in  order  to  reduce  the  tunneling  losses. 
The  sample  was  grown  by  molecular  beam  epitaxy  (MBE),  which  offers  high  accuracy 
for  the  design  parameters,  such  as  doping  and  layer  thickness,  since  we  have  disposition 

of  one  or  more  pure  materials  onto  the  crystal  wafer,  one  layer  of  atoms  at  the  time,  under 
ultra-high  vacuum.  The  preliminary  design  of  the  QWIP  predicted  an  IR  absorption  peak 
at  10.2  pm  [26]. 
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Figure  3.1:  Diagram  of  the  designed  and  simulated  QWIP  (From  [26].) 


C.  DESCRIPTION  OF  THE  COMPUTER  MODEL 

The  structure  used  in  simulation  is  a  InxGax_xAs  single  step  quantum  well  as 

shown  in  Figure  3.2. 
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Figure  3.2:  Potential  diagram  of  simulated InxGax_xAs  quantum  well  structure. 
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This  quantum  well  is  divided  into  four  regions.  The  percentage  of Inx  in  the  first, 

and  the  fourth  regions  is  zero,  while  for  regions  two  and  three  is  30  and  10  percent,  re¬ 
spectively.  Similarly  the  width  of  the  regions  one  and  four  is  considered  infinite,  while 
regions  two  and  three  are  40  Angstroms  wide  each. 

The  behavior  of  a  InxGax_xAs  single  quantum  well  was  investigated  in  the  three 

main  cases.  In  the  first  case  we  have  studied  the  InxGaXxAs  quantum  well  without  the 
application  of  an  electric  field,  in  the  second  case  the  electric  field  is  present,  and  finally 
in  the  third  case  the  electric  field  is  present  but  its  value  is  relatively  small,  which  needed 
a  special  treatment.  The  main  difference  between  these  three  cases,  is  that,  in  the  first 
case,  we  have  used  the  exponential  solution  of  the  Schrodinger  equation,  while  in  the 
case  of  large  and  small  bias  we  took  advantage  of  the  Airy  and  asymptotic  Airy  solutions, 
respectively  (Appendix  C). 

1.  Simulation  Parameters 

For  our  numerical  calculations,  the  values  of  the  effective  mass  m* ,  potential 
height  V  ,  and  energy  gap  E  for  the  InxGax_xAs  quantum  well,  were  employed  as: 


m  =  - 


0.028 m. 


-  + 


(0.067 m0) 


o  ) 


(3.1) 


where  m0  is  the  rest  mass  of  the  electron  and  %xj  is  the  mole  fraction  of  In  in  the  region 
j  .  The  bandgap  energy  Eg  (in  eV  )  of  the  four  layers  and  the  corresponding  barrier 
heights  F  can  be  obtained  by  the  following  formulas  according  to  [1]: 

Elg  =  EAg  =  1 ,519  +  1 .247  ( %Xj ) ,  (3.2) 

E\  =  E]  =  1 .5 19  - 1 . 102  ( %Xj ) ,  (3.3) 

Vj=0.62  (EJg-E;).  (3.4) 


Equations  (3.1)  through  (3.3)  characterize  the  physical,  structural  and  electronic 
properties  of  InxGax_xAs  under  specific  conditions. 
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D.  STEP  QUANTUM  WELL  WITHOUT  AN  APPLIED  ELECTRIC  FIELD 


1.  Bound  State  Energies 

The  bound  state  energies  are  calculated  using  the  transfer  matrix  method  (TMM), 
according  to  the  derivation  of  Subsection  E-2  of  the  Chapter  II.  The  energy  eigenvalues 
are  derived  by  the  condition m21  (£)  =  (),  given  in  Equation  (2. 14).  In  order  to  find  the 

values  of  E  that  satisfy  the  above  condition,  m12  (E)  was  evaluated  for  energies  from  the 
bottom  of  the  quantum  well  ( E  =  0  ),  to  the  top  ( E  =  Vl  =  V4  ),  using  small  increments  A E  . 
In  the  numerical  analysis,  the  m21  (  E )  derived  at  each  energy  was  multiplied  by  the  corre¬ 
sponding  value  at  the  previous  energy  and,  if  the  resulting  quantity  is  negative,  then 
m22  ( E )  had  crossed  a  zero.  This  procedure  allowed  us  to  detennine  the  bound  state  en¬ 
ergies  in  the  quantum  well. 


In  the  case  of  the  step  quantum  well,  there  are  four  layers,  which  require  three 
matrices  to  relate  the  coefficients  of  the  outermost  layers  as  given  in  Equations  (3.5)  to 
(3.7): 


ro 

rcA 

( 0  > 

f  ml. 

m  p  ^ 

A  > 

(m\1Dx'\ 

=  MX 
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Iaj 
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1^22  A  J 
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vA  y 

\m2\ 

m22; 
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vAy 


=  M3 


\D3J 


=m2m2m1 


AA  (m™  ^ 


vAy 


mu  mn 

321  321 

ym2\  m22  J 


vAy 


<2  A 
v<Ay 


(3.5) 

(3.6) 

(3.7) 


where  M j ,  j  =  1, 2,3  ,  is  given  by  Equation  (2.12)  and  niuv ,  m2Jv ,  nA , where  /i,v  =  1, 2  , 
are  the  matrix  elements  of  M]  ,the  resulting  elements  of  the  multiplication  of  M2M \ ,  and 
the  resulting  elements  of  the  multiplication  of  M3M2M \ ,  respectively. 


The  next  step  involves  normalization  of  the  bound  wavefunction,  which  allows 
the  detennination  of  Dx .  Normalization  implies: 


|  vPy.(x)vP*(x)  dx  =  1, 


—oo 


(3.8) 
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where  'F(x)  and  (V*(x)  are  the  wavefunction  and  its  complex  conjugate,  respectively. 

For  the  parameters  used  for  the  step  quantum  well,  there  was  only  one  bound  state 
in  the  well  as  shown  in  Figure  3.3,  alone  with  the  normalized  wave  function. 


X  (Angstroms) 

Figure  3.3:  InxGax_xAs  step  quantum  well  ground  state  energy  without  applied  bias. 

2.  Continuum  State  Energies 

The  responsivity  of  QWIP  detectors  is  strongly  dependent  on  the  location  of  the 
excited  state  [1,  17,  18].  To  avoid  the  suppression  of  photocurrent  due  to  tunneling,  the 
excited  state  is  usually  aligned  with  the  barrier  which  also  minimizes  the  thermionic 
emission  of  ground  state  electrons  [1].  It  is  possible  to  calculate  the  continuum  state 
wavefunctions  by  using  the  transfer  matrix  method,  with  proper  boundary  conditions.  If 
the  electron  is  incident  from  the  left  side  of  the  potential  well,  as  shown  in  Figure  3.2, 
then  the  coefficient  D4  =  0 ,  since  it  is  not  possible  to  have  a  left  moving  wave  in  the  right 
side  of  the  well. 

Using  the  transfer  matrix  method,  the  coefficients  of  the  fourth  and  the  first  re¬ 
gions  can  be  related  as: 


25 


=  M2M2Ml 


321  _ 321 


OTjf  m22 


The  coefficients  D,  and  C4  can  be  obtained  in  terms  of  the  matrix  elements  de¬ 


fined  earlier  as: 


D  —  —  _2i —  c 

1  _ 321  M  ’ 


(3.10) 


r  —  mlxC  +m321D  ~rnlxC  —  —  =  ^ 1  -±  171x2  ^21 

W  rnU  M  mll  M  321  321 


(3.11) 


Similarly,  the  transfer  matrix  method  gives  the  coefficients  of  the  continuum 
wavefunctions  for  the  second  and  third  regions: 


J  C  _  m^m2\ 

ai'-'i  ...321 


1 


(3.12) 


=MM, 


mu  mn 

m2l  m22 


21  ' 
m\  i  Q 


21  321 

2i  tn22wi2l 


(3.13) 


Figure  3.4  shows  the  dependence  of  coefficients  in  different  layers  as  a  function 
of  energy,  for  the  step  quantum  well  structure. 
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Figure  3.4:  Coefficients  of  the  continuous  states  vs.  energy  difference. 

Normalization  of  the  continuum  wavefunctions  cannot  be  executed  using  Equa¬ 
tion  (3.8),  since  these  wavefunctions  extend  to  infinity.  The  nonnalization  of  the  contin¬ 
uum  wavefunctions  are  usually  done  either  using  a  hypothetical  box  around  the  quantum 
well  [18],  or  using  momentum  space  5  function  normalization  as  described  in  [20].  Fig¬ 
ure  3.5  illustrates  100 continuum  wavefunctions  up  to  0.285  eV  above  the  ground  state 
using  an  arbitrary  scale.  The  red  curve  in  Figure  3.5,  illustrates  the  ground  state  wave- 
function.  The  oscillatory  behavior  of  continuum  wavefunctions  represents  the  free  mo¬ 
tion  of  in  the  continuum. 
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X-Angstroms 

Figure  3.5:  Continuum  wavefunctions  (blue)  and  ground  state  wavefunction  (red). 

E.  STEP  QUANTUM  WELL  STRUCTURE  UNDER  AN  APPLIED  ELEC¬ 
TRIC  FIELD 

1.  General 

As  described  in  Chapter  II,  it  is  necessary  to  use  Airy  functions  to  represent  wave- 
functions  when  the  quantum  well  is  under  bias.  The  implementation  of  Airy’s  function 
routines  in  a  computer  proved  to  be  a  relatively  complicated  task,  especially  when  the 
electric  field  is  small.  Thus,  a  “weak  field  limit”  was  placed  where  the  built-in  Matlab 
Airy  functions  give  numerical  overflows,  so  lower  biases  than  this  threshold  were  treated 
using  the  asymptotic  expressions  of  Airy  functions.  This  limit  was  found  to  be  about 
105  V/m  .  The  analytic  coding  used  for  this  case,  along  with  the  previously  examined 
cases  is  presented  in  the  Appendix  A. 

2.  Bound  State  Energies 

The  following  shows  the  calculation  of  bound  state  energies  in  the  step  quantum 
well  under  a  moderate  electric  field  of  106  V/m .  The  energy  eigenvalues  can  be  found  by 
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satisfying  the  condition  (2.52),  while  the  coefficients  of  Airy  functions  for  different  re¬ 
gions  can  be  obtained  as: 

C,  =  m\l\C2  =  =  ml,C4  =  0 ,  (3.14) 

D1  =  0 ,D2  =  mj3 , D3  =  m21,  D4  =  1,  (3.15) 

where  m3iv ,  m23v ,  in'23 ,  //,  v  =  1, 2  ,  are  elements  of  M3 ,  M2M3 ,  MjM2M3  matrices,  respec¬ 
tively  .  There  was  only  one  bound  state  in  the  well  as  shown  in  Figure  3.6: 


X  (Angstroms) 

Figure  3.6:  Ground  energy  and  wavefunction  in  the  step  quantum  well  under  an  elec¬ 
tric  field  of  106  V/m . 

In  the  case  of  negative  applied  electric  field  of  -106  V/m  the  ground  state  energy 
increases  compared  to  the  unbiased  quantum  well  as  shown  in  Figure  3.7. 
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Figure  3.7:  Ground  energy  and  wavefunction  in  a  InxGd\_xAs  quantum  well  under 

-106V/m  applied  field. 

It  is  well  known  that  the  bound  state  energies  in  a  quantum  well  do  not  shift  when 
measured  from  the  center  of  the  well  [1,  17].  Since  our  origin  is  at  the  edge  of  the  well, 
the  amount  of  the  energy  shift  can  be  estimated  using: 

u  fl 

shift  =  --L-L. —  s  (3.16) 

where  e  is  the  electron  charge,  F  is  the  applied  bias,  L  is  the  length  of  the  well  which 
confines  the  energy  state.  This  implies  that  for  a  106  V/m  applied  field,  the  shift  is  about 
-0.002  eV  from  the  unbiased  position,  and  for  -106  V/mthe  shift  is  about0.002  eV  . 

This  gives  an  energy  separation  of  0.004  eV  for  ±106  V/m  field  which  is  very  close  to  the 
simulated  value  of  0.005  eV  . 

In  the  case  of  application  of  an  electric  field  of  7  x  104  V/m ,  the  asymptotic 
expressions  of  Airy  functions  were  needed  to  obtain  the  ground  state  energy  and  the 
associated  wavefunction  are  illustrated  in  Figure  3.8. 
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Figure  3.8:  Ground  energy  and  wavefunction  in  a  quantum  well  under  7 x  104  V/m  . 


In  the  case  of  the  same  bias  applied  in  the  negative  direction,  the  program  gives 
that  the  ground  energy  Ex  lies  at  0.091  eV  ,  as  shown  in  Figure  3.9. 


X  (Angstroms) 


Figure  3.9:  Ground  energy  and  wavefunction  in  a  quantum  well  under -7  x  104  V/m . 
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Again,  the  difference  between  the  values  of  the  ground  energy  level  can  be  exp¬ 
lained  by  Equation  (3.18)  since,  for  an  7x  104  V/m  applied  field,  the  shift  is  about 
-7 xlO  4  eV  from  the  unbiased  position,  and  for  -7x  104  V/m  is  about 7  x  1CT4  eV .  This 
gives  an  energy  separation  of  0.00143  eV  for  ±7  x  104  V/m  field  which  is  very  close  to 
the  simulated  value  of  0.00137  eV .  Finally  the  lowest  bias  that  we  have  achieved  in  our 
simulation  was  ±7xl04  V/m  or  ±0.7  kV/cm  ,  which  is  much  smaller  than  typically  used 
in  QWIP  detector  operation  (10  kV/cm)  [1,  16]. 

3.  Continuum  State  Energies 


The  continuum  state  wavefunctions  are  obtained  by  assuming  the  electron  is  inci¬ 
dent  from  the  right  side  and  setting  the  following  boundary  conditions  for  the  coefficients 
of  the  outermost  layers: 

D,  =  0 .  (3.17) 


Then  applying  the  TMM,  the  coefficients  in  each  layer  can  be  obtained  as: 
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(3.22) 


For  the  normalization  of  the  continuous  wavefunctions,  the  technique  described  in 
the  Appendix  of  Ref.  [21]  was  used,  since  8  function  normalization  was  not  valid  due  to 
spatial  dependence  of  the  potential  energy. 
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The  absolute  values  of  the  coefficients  as  a  function  of  continuum  energy  are 
illustrated  in  Figure  3.10.  In  the  same  figure,  we  observe  the  oscillatory  behavior  of  the 
coefficients  due  to  the  resonances  formed  in  the  continuum  region,  due  to  constructive 
interference  of  incident  and  reflected  waves. 


x  104 


Energy 


Figure  3.10:  Absolute  value  of  coefficients  under  106  V/m  forward  bias. 

The  first  100  continuum  wavefunctions  of  the  tilted  well  are  illustrated  in  Figure 
3.11,  along  with  the  ground  state  wavefunction. 
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Figure  3.11:  Continuum  wavefunctions  (blue)  and  ground  state  wavefunction  (red)  for 
the  step  quantum  well  under  106  V/m  forward  bias. 


F.  BOUND-TO-CONTINUUM  TRANSITIONS  IN  A  QUANTUM  WELL 


1.  Fermi’s  Golden  Rule 

It  was  shown  by  Enrico  Fermi  [9],  that  the  transition  probability  IF  of  a  quantum 
system  from  an  initial  state  /  to  a  final  state  F  is  given  by: 
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w  =— y 


\V  hF 

F\v  p\  1  1 


5{Ef-Ei-Tico), 


(3.23) 


where  vPf,  and  lF/  are  the  final  and  the  initial  wavefunctions,  respectively,  EF  and  Ej 
the  associated  energies,  Vp  is  the  interaction  potential  and  fi  a>  the  interaction  photon 

energy.  The  photon  interaction  potential  in  dipole  approximation  is  given  by  [2]: 

f  T  A  V/2 

ZPe, 
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k2so£COC  J 


(3.24) 


where  e  is  the  electron  charge,  m  the  electron  effective  mass,  I  is  the  incident  photon 
flux,  £oands  characterize  the  electric  permittivity,  c  is  the  speed  of  light,  co  is  the 
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frequency  of  the  photon,  pe  defines  the  momentum  of  the  electron  and  s  is  the 

polarization  vector  of  the  photon. 

2.  Oscillator  Strength 

The  oscillator  strength  fosc  for  a  transition  from  an  initial  state  to  a  final  state  is 
defined  as  [1]: 

Lc=^-(EP-E,)(z)2  =^f(Ep-EI)\"'V](z)z'¥„(z)dz,  (3.25) 

n  n  Jzi 

where  EF  and  vFf  are  the  final  energy  and  wavefunction,  respectively,  El  and  'F/  are 

the  initial  energy  and  wavefunction,  m*  is  the  effective  mass,  ti  is  Plank’s  constant 
and  z  is  the  direction  of  growth. 

Generally  the  oscillator  strength  is  an  indication  of  how  strong  the  transition  is, 
which  directly  translates  to  the  performance  of  a  QWIP  detector. 

Figure  3.12  shows  the  calculated  oscillator  strength  for  the  step  quantum  well, 
under  zero  bias,  using  the  wavefunctions  derived  in  Subsection  D-2.  The  energy  was 
measured  from  the  ground  state. 
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Figure  3.12:  Oscillator  strength  versus  energy  for  bound-to-continuum  transitions  in  an 

unbiased  InxGax_xAs  step  quantum  well. 
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Under  an  applied  electric  field  the  oscillator  strength  can  be  calculated  using  the 
wavefunctions  derived  in  Subsection  D-3.  Figure  3.13  shows  the  oscillator  strength  for 
106  V/m  applied  bias. 


Figure  3.13:  Oscillator  strength  under  106  V/m  applied  bias  vs.  energy. 

The  oscillatory  behavior  of  the  oscillator  strength  is  due  to  the  constructive  inter¬ 
ference  of  incident  and  reflected  electron  waves,  from  the  tilted  left  quantum  well  barrier. 

This  oscillatory  behavior  becomes  denser,  as  the  strength  of  the  electric  field  was 
reduced,  as  illustrated  in  Figure  3.14.  As  the  bias  approaches  zero,  the  oscillations  be¬ 
come  closer  and  closer,  to  form  a  continuous  curve  as  shown  in  Figure  3.12. 
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Figure  3.14:  Oscillator  strength  under  7 x  104  V/m  applied  bias  vs.  energy. 

The  origin  of  oscillations  as  a  function  of  electric  field  can  be  understood  by  esti¬ 
mating  the  amount  of  energy  change  (AE )  required  for  obtaining  path  difference  between 

incident  and  reflected  electron  waves  ( Az )  by  one  wavelength.  This  can  be  estimated  by 
using  the  linear  potential  generated  by  the  electric  field  as: 

^  =  \e\F/±z.  (3.26) 

It  can  be  seen  from  Equation  (3.26)  that  as  the  electric  field  becomes  smaller,  the 
amount  of  energy  change  needed  to  get  one  wavelength  path  difference  reduces.  This 
implies  that  the  constructive  interference  occurs  at  a  rapid  rate,  as  the  energy  of  the  elec¬ 
tron  goes  up,  resulting  a  rapid  oscillation  of  the  oscillator  strength. 

3.  Absorption  Coefficient 

The  performance  of  photodetectors  can  be  conveniently  characterized  using  the 
absorption  coefficient  which  is  defined  as: 
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(3.34) 


a^hco) 


/z<yx(number  of  transitions  per  unit  volume  and  time) 
incident  energy  flux 


For  a  QWIP  detector  the  absorption  coefficient  can  be  expressed  using  the  density 
of  states  and  oscillator  strength  as  [1,2]: 
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(3.35) 


where  Nd  is  the  doping  density  in  the  well,  L  is  the  total  length  of  the  doped  region,  m*b 
and  Vb  are  the  barrier  effective  mass  and  potential  height,  respectively,  and  fi  co  is  the 
incident  photon  energy.  Figure  3.15  shows  the  oscillator  strength  and  the  density  of 
states  for  the  step  quantum  well  under  a  106  V/m  electric  field. 


Figure  3.15:  Oscillator  strength  and  density  of  states  vs.  energy. 

The  absorption  coefficient  was  calculated  for  the  step  quantum  well  assuming  the 
doping  density  to  be  1018  cm-3  and  the  result  is  shown  in  Figure  3.16.  The  calculated 
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absorption  coefficient  is  in  good  agreement  with  that  typically  observed  experimentally 
(□  800  cm-1)  [22]. 
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Figure  3.16:  Absorption  coefficient  vs.  photon  energy  under  106  V/m  applied  bias. 


G.  COMPARISON  WITH  EXPERIMENT 


The  step  quantum  well  structure  used  for  the  simulation  in  the  present  Chapter 
was  fabricated  and  characterized  by  Touse  [26]  and  Yeo  [27]  using  photocurrent  spec¬ 
troscopy.  In  this  section,  the  measured  responsivity  of  the  detector  is  compared  with  the 
simulated  responsivity  using  the  absorption  coefficient  and  outgoing  probability  current 
of  the  excited  electron. 

1.  Responsivity 

For  a  detector  the  nonnalized  voltage  or  current  signal  per  incident  optical  power 
is  called  responsivity  [24].  We  can  express  the  responsivity  R  as: 

R  =  rjj~g,  (3-36) 

hv 

where  77  is  the  quantum  efficiency,  e  is  the  electron  charge,  hv  is  the  average  photon 
energy  and  g  is  the  photoconductive  gain. 
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Taking  into  account  all  the  geometrical  and  optical  factors,  the  test  detector  re- 
sponsivity  can  be  estimated  using  [23]: 

I , .  fl  •  Idt  ■  R  f  ■ d-Sw 

R  —  — — —  — - - — ^ - ,  (3.37' 

P  V  ■  A  T  T 

1  det  v  ref  ^det  1  ZnSe  1  GaAs 

where  7det  refers  to  the  photocurrent,  Pdet  is  the  optical  power  incident  on  the  test  detector 
power,  Rref  is  the  responsivity  of  the  reference  detector,  d  ■  Sw  refers  to  the  reference  de¬ 
tector  surface  area,  Vref  is  the  generated  by  the  reference  detector  voltage,  Adet  is  our  own 
detector  surface  area,  TZnSe  and  TGaAs  are  the  transmittances  of  Zinc  Selenide  and  Gallium 
Arsenide,  respectively,  of  the  cold  head  where  the  detector  was  situated. 

During  the  experiment  the  temperature  was  kept  fixed  at  40K,  so  the  photocur¬ 
rents  for  the  test  detector  under  different  applied  electric  fields  were  measured,  and  the 
responsivities  calculated  using  Equation  (3.37). 

The  experimental  responsivity  as  a  function  of  wavelength,  for  the  step  quantum 
well  detector,  for  a  set  of  biases  is  illustrated  in  Figure  3.17.  The  shifting  of  the  peak  po¬ 
sition  in  Figure  3.17  is  mainly  due  to  the  linear  Stark  effect  in  step  quantum  well  as  ex¬ 
plained  in  Chapter  II. 


wavelength  nm 

Figure  3.17:  Experimental  responsivity  vs.  wavelength  at T  =  40  K  (From  [27].) 
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2.  Absorption 

The  normalized  absorption  coefficients  calculated  using  the  approach  described  in 
Subsection  F-3  of  the  present  Chapter,  as  a  function  of  wavelength  for  the  bias  voltages 
used  in  the  experiment  are  shown  in  Figure  3.18.  The  shifting  of  the  peak  position  can  be 
clearly  seen  from  Figure  3.18. 


Figure  3.18:  Computer  simulation  of  absorption  vs.  wavelength  for  the  step  quantum 

well  structure. 


3.  Photocurrent 

The  photocurrent  density  can  be  estimated  by  taking  the  product  of  quantum  effi¬ 
ciency  with  the  outgoing  current  density  of  the  excited  electron.  The  current  density  as¬ 
sociated  with  a  quantum  state  can  be  calculated  using  the  wavefunction  [24],  as: 


J  = 


Ti\e\ 


2  mi 


¥ 


(TV  (TV 


dx 


dx 


(3.38) 


where  ti  is  the  reduced  Plank’s  constant,  m  is  the  electron  mass,  x  is  the  direction 
where  the  current  is  evaluated  and  'F  is  the  wavefunction. 

Assuming  the  wavefunction  has  the  following  fonn: 
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CAi  +DBi, 


(3.39) 


and  after  some  algebraic  effort,  the  quantity  in  the  bracket  of  Equation  (3.38)  can  be  ob¬ 
tained  as: 


qr  _  xp  =  (c*D- D*C) ( AiBi  -  AiBi )  , 


dx 


(3.40) 


Using  the  Wronskian  (W)  of  the  two  Airy  functions  [40]  the  second  bracket  in 
Equation  (3.40)  can  be  further  simplified  as: 


W  = 


Ai  Bi 
Ai  Bi 


=  AiBi  -  Ai  Bi  =  — . 

n 


(3.41) 


The  current  density  given  by  Equation  (3.38)  can  be  expressed  using  the  result  in 
Equation  (3.41)  as: 

h\e\ 


J  = 


-  (C*D-D*C), 


(3.43) 


2  ni  m\p\ 

where  fd  is  given  by  Equation  (2.20),  and  appears  in  the  denominator  of  Equation  (3.43) 
due  to  the  argument  of  the  Airy  functions  of  Equation  (3. 39), which  according  to  Equation 
(2.39)  is: 
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(3.44) 


Equation  (3.43)  represents  the  total  current  density  which  should  be  zero,  since 
the  probability  density  is  time  independent  [20].  This  requires  that  the  coefficients 
C  and  D  must  be  related  as: 


C*  xD  =  D*  xC  .  (3.45) 

In  our  simulations,  it  was  found  that  this  condition  holds  for  all  the  regions  of  the 
quantum  well  structure.  This  indicates  that  the  current  density  given  in  Equation  (3.43) 
consists  of  two  identical  components,  (a)  outgoing  current  and  (b)  incoming  current, 
which  cancel  to  give  zero  net  current. 
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In  order  to  estimate  the  photocurrent  after  excitation  of  the  electron  from  the 
ground  state,  it  is  necessary  to  extract  the  outgoing  component  from  the  current  density 
that  is  given  in  Equation  (3.43).  The  incoming  and  outgoing  components  of  the  current 
are  not  obvious  in  the  case  of  the  linear  combinations  of  Airy  functions,  while  they  are 
usually  obvious  for  exponential  wavefunctions  obtained  for  the  zero  bias  case. 

However,  it  can  be  shown  that  the  linear  combinations  shown  in  Equations  (3.46) 
and  (3.47)  produce  pure  outgoing  and  incoming  currents  by  substituting  in  Equation 
(3.38): 

^out=Ai-iBi,  (3.46) 

Tf  =  Ai  +  i  Bi .  (3.47) 


The  magnitudes  of  current  density  corresponding  to  the  outgoing  and  incoming 
wavefunctions  are  given  by: 
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(3.48) 

(3.49) 


Under  a  positive  bias,  the  photoexcited  electrons  moves  to  the  right  of  the  quan¬ 
tum  well  and  the  wavefunction  in  outermost  right  region,  given  in  Equation  (3.50)  can  be 
used  to  determine  the  outgoing  current. 

'F,  =C4Ai  +  D4Bi,  (3.50) 

where  C4  and  Z)4  are  the  coefficients  defined  in  Subsection  D-2  of  Chapter  II . 


Equation  (3.50)  can  be  rewritten  using  the  outgoing  and  incoming  wavefunctions 
given  in  Equations  (3.46)  and  (3.47)  as: 


'F  =  C'F  +  D'F 

where  the  constants  C  and  D  are  related  to  the  coefficients  C4  and  D4 , since: 

r,  _C4-iDA 
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The  corresponding  currents  are  given  by: 
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(3.54) 

(3.55) 


The  photocurrent  density  was  calculated  using  the  outgoing  current  given  in 
Equation  (3.54)  and  the  quantum  efficiency  rj  =  o.L  as: 


h 

C4  -  i  D4 

n  m  \P\ 

2 

(3.56) 


where  a  is  the  absorption  coefficient  and  L  is  the  total  length  of  the  quantum  well  struc¬ 
ture. 


Figure  3.19  shows  the  calculated  photocurrent  for  the  biases  used  in  the  experi¬ 
mental  measurement  of  the  step  quantum  well  detector. 


F igure  3.19:  Calculated  photocurrent  vs .  wavelength. 
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Figure  3.20  compares  the  absorption  and  normalized  photocurrent  density  with 
the  measured  normalized  responsivity  at  0.81  V  forward  bias.  It  can  be  observed  from 
this  figure  that  the  simulated  photocurrent  density  provides  a  good  description  of  the  ex¬ 
perimental  observations.  The  absorption  is  found  to  be  higher  than  the  photocurrent  den¬ 
sity  in  the  longer  wavelength  regime.  This  is  mainly  due  to  the  difficulty  of  extraction  of 
electrons  after  excitation  due  to  the  triangular  barrier  on  the  right  side  of  the  biased  quan¬ 
tum  well. 


Figure  3.20:  Measured  responsivity,  simulation  photocurrent  density  and  absorption  vs. 

photon  wavelength  for  0.81  V  forward  bias. 

H.  SUMMARY 

In  this  chapter,  we  presented  the  calculation  of  absorption  and  photocurrent  den¬ 
sity  for  the  step  quantum  well  detector  under  an  external  bias.  The  calculated  results 
were  compared  with  the  experimental  data  and  found  to  provide  a  good  agreement  which 
validated  the  accuracy  of  the  model  employed.  In  the  following  chapter,  we  will  present 
a  simulation  code  based  in  EKF  for  angle  only  tracking  using  an  IR  sensor.  Ideas  such  as 
initialization  using  two  synchronous  bearing  and  EKF  angle  tracking  will  be  developed 
and  finally,  selected  results  will  be  presented. 
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IV.  PERFORMANCE  OF  QUANTUM  WELL  INFRARED 
PHOTODETECTORS  IN  TARGET  TRACKING 


A.  INTRODUCTION 

Quantum  well  infrared  photodetectors  constitute  a  tool  for  target  tracking  since 
they  offer  the  advantage  of  passivity.  As  we  have  discussed  in  the  previous  chapters,  tak¬ 
ing  advantage  of  our  ability  to  design  multi-wavelength  tunable  infrared  detectors,  we 
can  combine  them  with  batch  processing  or  recursive  type  algorithms  in  order  to  track 
targets  passively. 

In  this  chapter  we  simulate  in  Matlab  the  scenario  of  bearing-only  tracking  of  a 
target  using  a  recursive  algorithm  based  on  the  Kalman  filter.  We  investigate  the  bearing- 
only  tracking  problem  in  two  main  cases,  namely  tracking  a  non-maneuvering  target  and 
tracking  a  maneuvering  target  using  either  one  or  two  sensors.  Since  we  study  only  the 
tracking  problem,  we  suppose  that  the  detectors  are  tuned  in  such  way  that  they  provide 
the  highest  contrast  and  avoid  the  extensive  clutter  so  they  generate  relatively  reliable 
bearings.  More  analytic  information  about  the  coding  used  can  be  found  in  Appendix  B. 
Finally  the  BOT  problem  is  considered  trivial  since  limited  information  is  available 

B.  DESCRIPTION  OF  THE  MODEL 

The  simulation  model  and  the  encounter  geometry  are  illustrated  in  Figure  4. 1 . 
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Sensor  1  is  located  at  position  (  0,0)  of  the  local  coordinate  system,  while  sensor 

2  was  placed  at(10km,0) .  These  sensors,  since  they  are  IR,  can  make  measurements  of 

the  bearing  of  the  target  only.  The  first  step  of  the  bearing-only  tracking  using  an  ex¬ 
tended  Kalman  filter  is  the  initialization  of  the  filter.  Initialization  in  the  case  of  EKF  is 
sufficient  if  we  estimate  two  positions  of  the  target.  This  is  possible  either  by  taking  ad¬ 
vantage  of  the  observer  and  target  relative  trajectory  in  the  case  that  we  use  one  sensor,  or 
by  assuming  synchronous  multiplatform  bearings.  BOT  using  asynchronous  bearings  will 
not  be  examined.  So  without  the  loss  of  generality  in  our  case,  we  assume  two  stationary 
IR  observers  and  synchronous  initialization  and  tracking  bearings. 

C.  INITIAL  ESTIMATE  OF  THE  TARGET  POSITION 

The  initial  estimate  of  the  target  position  is  determined  using  two  sets  of  time  syn¬ 
chronous  bearings.  We  have  used  the  fact  that,  in  order  to  initialize  our  filter,  we  do  not 
need  to  have  the  exact  position  of  the  target,  only  a  reliable  estimate  of  its  position.  The 
simplified  geometrical  problem  is  illustrated  in  Figure  4.2. 


We  suppose  that  the  first  set  of  the  synchronous  bearings  from  our  IR  sensors  are 
6X  and  02  as  shown  in  Figure  4.2.  The  reliability  of  the  bearings  depends  on  the  quality  of 
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the  sensors,  the  sensors-target  geometry  and  the  distance  between  the  sensors  and  the  tar¬ 
get.  Since  we  examine  the  general  case  of  the  problem  of  initialization,  we  suppose  that 
the  sensors  we  use  give  a  random  error  in  bearing  of  maximum  magnitude  of  approxi¬ 
mately  2°  .  This  error  is  considered  very  high  for  sensors  designed  to  provide  only  the 
bearing,  and  it  can  be  positive  or  negative. 

So,  examining  sensor  1,  the  angular  real  position  of  the  target  must  be  between: 

#1 ~^(Xtn,e>ytrue)^0l>  C4’1) 

where 

#f  =  0l-  |rand(0, 1)  •  2°| ,  (4.2) 

9+  =01+|rand(O,l)-2°|,  (4.3) 

while  bearings  9~  and  9]  are  shown  in  Figure  4.2  with  the  blue  dashed  line. 


Treating  the  sensor  2  case,  using  the  methodology  used  for  sensor  1,  we  observe 
that  the  target  must  be  inside  the  quadrilateral  B CDK  .  In  order  to  find  the  coordinates  of 
BCDK  we  apply  Euclidian  geometry  to  each  triangle  of  Figure  4.2.  So,  examining  trian¬ 
gle  OB  A  and  applying  the  sine  law,  we  get: 

sin(0,+)  =  sin(;r-fl2+)  =  sin(#j)  ^ 

AB  OB  OA  ’  1  ‘  ' 


where  AB ,  OB ,  OA  are  the  lengths  of  the  sides  of  the  triangle  as  shown  in  Figure  4.2 
while  9\  refers  to  the  third  angle  of  the  triangle,  namely  angle  OB  A  .  From  Equation 
(4.4)  we  can  find  the  length  OB  as: 


OB  =  OA 


sin  (n-6t) 


sin  (#f) 

where  OA  is  the  distance  of  the  two  sensors  while  angle#'  is  obtained  since: 


(4.5) 


0\  =n-0+  -  {n-  02 )  =  -  #+ . 


(4.6) 


Obtaining  the  length  of  OB  we  can  find  the  coordinates  using: 

sin(#2+)cos(#1+) 


x,  =  OB  cos(#j+)  =  OA  - 


sin(#2+  -#,+  ) 


(4.7) 
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yi=oasin(g>)  =  Qjsin(^)sin(^) 
1  1  sin(6)2+  -6+) 


(4.8) 


We  observe  that  finding  the  coordinates  of  point  B  involves  only  the  sensor  dis¬ 
tance  and  the  uncertainty  bearings  obtained  by  the  two  sensors.  Following  the  same  pro¬ 
cedure  the  coordinates  of  point  C  are  given  by: 


x2 


=  OCcos(d^)  =  OA 


sin(6,y)cos(6l1+) 

sin(^2"-^,+) 


(4.9) 


y2  =  OB  sin(6l1+)  =  OA 


sin(6h)sin(#j+) 
sin (6*2  -  Of) 


(4.10) 


For  the  coordinates  of  point  D  our  calculations  give: 


=  ODcos(Of)  =  OA 


sin(#r)cos(#f ) 
s  i  n  ( 6*2  -  0X  ) 


while  for  point  K 


y3  =  OD  sin(6*f  )  =  OA 


sin(02- )sin(0f) 
sin  (Of -Of) 


x4 


=  OK  cos (Of)  =  OA 


sin(#2+)cos(#f) 
sin(6l2+  -Of) 


y4  =  OK  sin(#f  )  =  OA 


sin(6,2+)sin(6*f  ) 
sin(#2+  -Of) 


(4.11) 

(4.12) 


(4.13) 

(4.14) 


The  next  step  involves  determination  of  the  maximum  between  the  two  diagonal 
lengths  BD  and  Of .  This  determination  is  used  for  initialization  of  the  covariance  ma- 


Figure  4.3:  Encounter  geometry  for  the  determination  of  the  diagonal  lengths. 
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The  diagonal  length  can  be  extracted  either  using  the  angle  (f)  of  the  orthogonal 


triangle  MCK  whose  inverse  tangent  is: 


<f)  =  tan-1 


x4  -x2 


(4.15) 


yyi-y*) 

while  the  length  CK  according  to  this  approach: 

|(*4-*2)| 


CK  = 


(4.16) 


sin 


tan- 


Vyi-y*j 


or  using  the  well  known  approach  for  finding  length: 


BD  =  ^](xl-x3)  +(yl-y3)  .  (4.17) 

Since  the  target  can  be  anywhere  inside  the  quadrilateral  we  assume  without  the 
loss  of  generality  that  it  is  located  at: 

-  V  -|-  V  -|-  v 

(4.18) 


x,  = 


yt  = 


x3  +  x2  +  x3  +  x4 


Ti  + v2+y3+y4 


(4.19) 


We  have  implemented  the  above  considerations  in  Matlab  and  we  have  produced 
Figure  4.4,  for  the  target  initialization: 


meters 


Figure  4.4:  Matlab  initialization  using  two  sets  of  synchronous  bearings. 
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The  position  of  the  target  during  the  first  and  the  second  time  that  the  bearings 
took  place  is  illustrated  in  Figure  4.5: 
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Figure  4.5:  The  position  of  the  target  is  located  inside  the  uncertainty  bearing  box. 

D.  EXTENDED  KALMAN  FILTERING  CONSIDERATIONS 
1.  General 

The  problem  considered  in  the  bearing-only  tracking  is  that  the  target  motion  is 
modeled  in  linear  Cartesian  coordinates,  while  the  measurements  are  in  polar  coordinates. 
We  begin  by  assuming  the  following  non-linear  discrete-time  system: 

x(k  +  \)  =  f(k,x{k))  +  g(k,x(k))vk, 
z(k)  =  h{k,x(k,y)  +  cok. 


where  x  provides  the  target  position  data,  z  is  the  measurement,  /  and  h  are  non-linear 
vector-valued  functions,  g  is  a  non-linear  matrix  valued  function  and  vk ,  cok  are  uncor¬ 
related  Gaussian  processes. 

According  to  [28-30],  it  can  be  proved  that  the  prediction  using  an  extended  Kal¬ 
man  filter  is  given  by: 

**+i|*=/(M*|Jt)>  (4-2l) 

Pk+V=FkW  +  GkQkG'k,  (4.22) 
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where  xk+V[k  is  the  state  vector  at  time  k  + 1  given  its  value  at  time  k ,  Pk+1,k  is  the  cova¬ 
riance  prediction  matrix  at  time  k  + 1  given  its  value  at  time  k  ,  Fk  and  Fk  are  the  gradient 
and  the  transpose  gradient  of  f(  k,xkk  ]j ,  Gk  and  G'k  are  the  zero  order  term  for  g{k,xk]k ) 
evaluated  at  xklk  and  its  transpose,  respectively  and,  finally,  Qk  is  a  matrix  related  to  the 
noise  inserted  in  the  system. 

We  define  Q,.  according  to  [27]  as: 

Qk  =  q  '-Q,  (4.23) 

where  Qt  is  the  sampling  interval  matrix  and  q  arbitrary  chosen  parameter  of  the  design 
[30]. 

The  update  of  the  measurement  according  to  [28,  29],  can  be  expressed  as: 

•Tt+i|jt-n  =  ^k+\\k  +  ^k+i  (za-+ i  —  h  {jc  + 1,  xk+V[k )  j ,  (4.24) 

P^i  =  (I- ~ Kk+lHk+l ) Pk+Mk  (/ - Kk+1Hk+l )'  +  Kk+1Hk+lK'k+1 ,  (4.25) 

where  zk+l  is  the  measurement  at  time  k  +  \  ,1  stands  for  the  associated  to  the  problem 
identity  matrix,  Hk+l  is  the  gradient  of  h[k  +  \*xk  U  j  and  Kk+l  is  defined  as: 

Kk+ 1  =  Pk+\\kH'k+\  {jH k+\Pk+\\kH'k+\  +  ^k+l  )  ’  (4-26) 

where  Rk+1  is  associated  to  the  uncertainty  of  the  measurement  at  time  k  + 1 . 

2.  Bearing  Tracking  Expansion  of  EKF 

The  initialization  in  the  EKF  BOT  case  provides  the  vector  of  the  two  positions  of 
the  target  after  the  second  set  of  the  synchronous  bearings.  This  vector  is: 

x  =  (x2y2xlyl)T,  (4.27) 

where  x2 ,  y2 ,  x, ,  y,  are  the  estimated  possible  positions  of  the  target  obtained  during  the 

initialization  process  by  the  two  sets  of  the  tautochronous  bearings.  Since  our  state  vector 
is  in  the  form: 

x  =  (xvxy  vy)T ,  (4.28) 
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where  x,y  is  the  position  of  the  target  and  vx ,  vy  the  associated  velocity,  we  use  matrix  D 
to  obtain: 

x  =  Dx,  (4.29) 

where  the  matrix  D  is  defined  as: 

10  0  0 

l/dt  0  -1/dt  0 

0  10  0 

0  l/dt  0  -l/dt 

where  dt  is  the  time  between  the  measurements. 

The  initialization  of  the  covariance  matrix  P  involves  the  maximum  diagonals  of 
the  quadrilaterals  and  it  is  given  according  to  [29,  30]  by: 

d2  0  0  0 

0  d\  0  0 

0  0  d2  0 

0  0  0  dl 


(4.31) 


where  dx  is  the  maximum  distance  of  the  diagonals  of  the  quadrilateral  fonned  by  the 
first  set  of  bearings,  while  d 2  refers  to  the  second.  During  the  EKF  update  and  for  the 
rest  of  the  calculation  the  covariance  matrix  P  is  obtained  by  Equation  (4.25). 

The  bearing-only  measurements  are  9  angle  measurements,  so  they  can  be  ex¬ 
pressed  by: 

z(k^  =  9{x{k)^  +  cok.  (4.32) 

where 


#(x(&))  =  arctan 


(4.33) 


The  gradient  of  h(x)  =  9(x)  for  the  bearing-only  measurements  can  be  expressed 


by  the  following  matrix: 


d6(x ) 

dx 


y 

x2  +  y2 


x2  +  y2 


(4.34) 


where  xand  y  are  obtained  by  Equation  (4.27). 
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3.  Maneuvering  Target  Considerations 

According  to  [29]  and  [30],  in  the  maneuvering  case  target,  the  equations  of 
measurement  and  update  do  not  change.  However,  the  transition  matrices  should  be  cal¬ 
culated  as: 


pu)  = 


sin(ruA) 

co 

cos(ruA) 

l-cos(ruA) 

co 

sin(ruA) 


0 

0 

1 

0 


l-cos(ryA) 

co 

-sin(<uA) 

sin(<z>A) 

co 

cos(ruA)  j 


(4.35) 


where  co  is  the  turn  rate  and  A  is  the  sampling  interval.  The  detennination  of  F(  j)  is 
possible  by  application  of  simple  target  motion  analysis  for  a  rotating  mass  around  a 
complete  circle.  The  turn  can  be  either  clockwise  or  anticlockwise. 

4.  Trajectory  Error  and  Measurement-Angle  Chi-square  Formulation 
The  error  is  detennined  by  the  difference  between  the  trajectory  that  the  EKF  pre¬ 
dicts  and  the  real  trajectory  of  the  target: 

vA  f  r  A 

(4.36) 


(x  3 

error 

(x  \ 
Areal 

\  y error  ) 

\y real ) 

where  x 


error  5  T error 


are  the  difference  between  the  predicted  by  the  EKF  x ,  y  and  the  real 


coordinates  xreal ,  yreal  of  the  position  of  the  target.  As  a  consequence  the  total  distance 


error  is  obtained  by: 


Total  Error  = 


V  y error  ) 


(  ^ error  T error  )  * 


(4.37) 


The  Chi-Square  values  of  the  angle  allow  us  to  observe  the  angle  variations  dur¬ 
ing  the  measurements.  These  values  according  to  [27]  are  obtained  by: 

chi2  =  (z-z)'  ■(H-P-H'  +  Ry1-(z-z),  (4.38) 

where  z  is  defined  by  Equations  (4.32)  and  (4.33),  z  is  the  non-linear  angular  measure¬ 
ment,  H  is  the  gradient  defined  by  Equation  (4.34),  P  is  the  covariance  matrix  and  R 
involves  the  error  of  the  measurement.  In  every  sampling  interval  the  chi2  value  changes 
since  the  measured  angle  of  the  target  is  different. 
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E. 


PROBLEM  FORMULATION  AND  RESULTS 


The  mathematical  considerations  above  are  the  basis  of  tracking  a  moving  target 
both  in  a  straight  line  and  during  a  turn  using  either  one  or  two  sensors,  with  or  without 
the  presence  of  measurement  noise.  For  the  purpose  of  this  thesis  we  simulated  each  case 
separately  in  Matlab,  although  only  selected  cases  will  be  presented.  Analytic  informa¬ 
tion  for  the  all  the  cases  can  be  found  in  Appendix  C.  In  our  simulation  we  assumed  that 
the  position  of  the  sensors  and  the  bearing  measurements  were  according  to  Figure  4.2. 
The  initialization  always  takes  place  using  two  sensors  since  these  are  assumed  station¬ 
ary,  but  it  is  possible  to  use  only  one  if  a  sensor  motion  is  assumed.  In  this  case  the 
mathematical  equations  should  be  corrected  in  such  way  that  the  relative  trajectory  of  the 
target-observer  is  taken  under  consideration  [29]. 

1.  Non-Maneuvering  Target  BOT  using  One  or  Two  IR  Sensors  in  Noise 
Presence 

In  the  case  of  a  non-maneuvering  target  the  initialization  of  the  filter  takes  place 
using  both  sensors,  while  during  the  tracking  we  can  use  either  one  or  two  without  any 
difference  in  the  results.  The  geometry  is  illustrated  in  Figure  4.6. 
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Figure  4.6:  True  position  and  track  position  of  straight  line  moving  target  using  two 
IR  sensors  in  the  presence  of  measurement  noise. 
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The  error  during  the  tracking  and  the  values  of  the  measurement  angle  are  pre¬ 
sented  in  Figure  4.7. 


Sample  number 

Chi-Squared  Values  for  Measurement  Association 


Figure  4.7:  Track  position  error  and  chi-squared  values  vs.  sample  number  for  non¬ 
maneuvering  target  -  two  IR  sensors  in  the  presence  of  measurement  noise. 


In  Figure  4.6  we  observe  that,  when  we  use  both  sensors  for  initialization  and 
tracking,  the  filter  tracks  the  flying  target  almost  from  the  starting  point  of  our  simulation 
run  until  the  final  point.  We  have  selected  a  geometry  for  the  simulation  which  does  not 
offer  advantage  to  the  EKF  since  the  target  flies  away  from  the  sensors  and  the  bearing 
angles  are  not  perpendicular.  Despite  these  facts,  the  filter  error  generally  remains  under 
200  mfor  most  of  the  simulation  time,  although  there  is  an  interval  of  40  samples 
where  the  error  increases. 

In  Figure  4.7  we  also  observe  that  the  values  of  the  chi-squared  variable  remain 
under  0-20  which  means  that  the  effort  of  the  filter  to  correct  the  bearing  angle  and, 
hence,  to  obtain  a  more  accurate  position  of  the  target  is  minimal.  This  means  that  the 
filter  in  unaware  that  there  is  a  difference  between  the  tracking  position  and  the  real  posi¬ 
tion  of  the  target.  This  result  was  expected  since  we  are  simulating  a  filter  that  tracks  only 
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with  angle  measurements,  so  in  a  case  when  the  EKF  is  used  with  radar  sensor  that  pro¬ 
vides  range  information,  these  effects  are  eliminated. 

2.  BOT  of  Low  Maneuvering  Target  using  IR  Sensors  in  Noise  Presence 
In  this  case  the  target  turn  was  simulated  according  to  Equation  (4.35).  We  will 
examine  two  cases  in  this  category.  The  first  case  involves  a  smooth  target  turn  of  not 
more  than  3  g ,  while  in  the  other  case  the  target  will  execute  a  sharp  high  g  turn.  The  ge¬ 
ometry  for  the  filter  and  the  target  for  the  first  case  are  shown  in  Figure  4.8. 


x  -|q4ES0  Tracking  single  Sensor, Noise  presence, q2  =  1000,  accel  =3g 


Figure  4.8:  True  position  and  track  position  of  a  low  g  turning  target  using  one  IR 

sensor  in  the  presence  of  measurement  noise. 

From  Figures  4.7  and  4.8  above  we  conclude  that  as  long  as  the  target  flies  on  a 
straight  course  the  single  sensor  filter  responds  positively.  When  the  target  starts  its  turn 
we  observe  a  slow  response  time  from  the  EKF  which  results  in  the  loss  of  the  target  for 

almost When  the  filter  realizes  that  a  turn  has  occurred,  it  turns  sharply  seeking  for 
the  true  target  track.  The  track  subsequently  generated  by  the  EKF  parallels  the  true  tar¬ 
get  path,  matching  the  true  bearing  closely  but  missing  in  range.  So  the  filter  tracks  the 
target  path  only  in  bearing.  This  behavior  appears  because  we  applied  the  EKF  algorithm 
to  track  a  turning  target  using  only  one  sensor  so,  once  the  sensor  loses  the  real  position 
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of  a  target  can  only  recover  the  target’s  bearing,  but  not  the  distance.  The  error  and  the 
chi-squared  graphs  presented  in  Figure  4.9  supports  the  previous  arguments. 
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Figure  4.9:  Track  position  error  and  chi-squared  values  vs.  sample  number  for  low  g 

maneuvering  target  -  single  IR  sensor  in  the  presence  of  measurement  noise. 

In  Figure  4.9  we  observe  that  using  both  sensors  for  initialization  and  a  single 
sensor  for  tracking  results  in  an  increasing  error  to  almost  700  meters  before  the  target 
starts  its  turn.  The  filter  tracks  the  target  almost  from  the  starting  point  of  our  simulation 
run  until  the  final  point.  The  tracking  in  this  case  is  characterized  as  low  quality  since  the 
true  position  of  the  target  after  the  turn  is  never  obtained.  At  sample  28  ,  as  we  clearly  ob¬ 
serve,  the  error  starts  growing  and,  at  the  end  of  the  simulation,  it  becomes  almost 
30  km  .We  also  observe  that  the  values  of  the  chi-squared  remain  under  0.10  when  the 
target  is  out  of  the  turn  loop  and  close  to  0.50  radians  when  the  target  is  inside  the  turn. 
The  fact  that  the  values  are  kept  low  means  that  the  filter  assumes  that  it  is  tracking  the 
target  accurately  after  the  turn.  This  is  true  in  bearing  but  not  in  range.  Again  since  we 
combine  the  EKF  with  a  single  angle  tracking  sensor,  the  algorithm  is  unable  to  regain 
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the  real  position.  The  tracking  results  become  worse  when  the  simulation  involves  higher 
random  errors  and  sharper  turns,  as  we  will  present. 

To  improve  the  performance  in  the  case  of  a3g  turn,  we  employed  both  sensors  in 
the  tracking  effort.  The  results  are  shown  in  Figures  4. 10  and  4. 1 1 . 


x  -|q4  BO  Tracking  two  Sensors, Noise  presence, q2  =  1000,  accel  =3g 
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Figure  4.10:  True  position  and  track  position  of  a  low  g  turning  target  using  two  IR 
sensors  in  the  presence  of  measurement  noise. 
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Track  Position  Errors:  BO  Tracking  two  Sensors, Noise  presence, q2  =  1000,  accel  =3g 


Figure  4.1 1 :  Track  position  error  and  chi-squared  values  vs.  sample  number  for  low  g 
maneuvering  target  -  two  IR  sensors  in  the  presence  of  measurement  noise. 

We  observe  that  using  the  second  sensor  has  improved  the  tracking  conditions; 
namely,  it  has  reduced  the  resulting  error  from  30  km  to  almost  4  km .  If  the  simulation 
was  free  to  run  for  a  longer  sampling  time,  the  real  position  and  the  tracking  position 
would  eventually  converge.  At  the  same  time  the  chi  squared  values  were  found  almost 
the  same. 

3.  BOT  of  High  Maneuvering  Target  using  IR  Sensors  in  Noise  Pres¬ 

ence. 

We  present  two  sub-cases  in  BOT  of  a  high  maneuvering  target.  Firstly,  the  concept 
of  EKF  angle  tracking  using  one  sensor.  In  this  case,  a  turn  of  no  more  than  9  g  is  pre¬ 
sented  since,  with  the  particular  geometry  and  with  the  particular  initialization  bearings 
and  sampling  time,  higher  turns  result  in  a  breakdown  of  the  EKF  algorithm.  So  employ¬ 
ing  a  9  g  turn  of  a  target,  which  is  moving  with  velocity  339  m/s ,  results  to  the  geometry 
of  Figure  4.12. 
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x  -|o4!30  Tracking  single  Sensor, Noise  presence, q2  =  1000,  accel  =9g 
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Figure  4.12:  True  position  and  track  position  of  a  high  g  turning  target  using  one  IR 
sensor  in  the  presence  of  measurement  noise. 


We  observe  that  the  single  sensor  combined  with  the  EKF  tracks  the  target  relia¬ 
bly  when  we  have  a  straight-line  motion,  but  when  the  turn  takes  place  it  follows  with  a 
time  delay  and  finally  manages  to  recover  the  correct  bearing  tracking.  The  associated 
error  and  the  values  of  chi-squared  are  illustrated  in  Figure  4.13. 
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Figure  4.13:  Track  position  error  and  chi-squared  values  vs.  sample  number  for 
high  g  maneuvering  target  -  single  IR  sensor  in  the  presence  of  measurement  noise. 

Figure  4.13  shows  that  the  error  is  kept  small  before  the  target  starts  the  turn,  it 
grows  to  almost  8  km  during  the  turn,  then  it  reduces  due  to  the  convergence  of  the  real 
and  the  estimate  position  but  finally  it  grows  due  to  the  divergence  of  the  two  tracks.  The 
chi-squared  values  are  also  higher  compared  to  the  other  cases  due  to  the  high  g  turn,  so 
we  observe  that  the  values  are  over  1  during  the  turn,  while  for  the  half  of  the  simulation 
time,  they  are  close  to  0. 1 . 

Second,  in  order  to  obtain  more  reliable  results  we  employ  the  second  sensor  in 
the  tracking  effort.  We  expect  the  performance  of  our  system  to  be  improved  similar  to 
the  results  of  the  low  g  turning  target.  In  fact,  the  tracking  results  turned  out  to  be  better 
than  expected,  meaning  that  the  IR  sensors  combined  with  the  EKF  algorithm  track  the 
target  with  small  error.  The  explanation  of  the  reliable  tracking  is  that  in  this  case  the  tar¬ 
get  follows  a  reverse  course  towards  our  system.  This  means  that  the  distance  to  our  sen¬ 
sor  system  is  reduced,  so  the  suggested  rule  that  the  tracking  distance  should  not  be  much 
higher  than  the  sensor  distance  is  satisfied  [30].  The  results  and  the  encounter  geometry 
are  shown  in  Figures  4.14  and  4.15. 
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Figure  4.14:  True  position  and  track  position  of  a  high  g  turning  target  using  two  IR 
sensors  in  the  presence  of  measurement  noise. 
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Figure  4. 15:  Track  position  error  and  chi-squared  values  vs.  sample  number  for  high  g 
maneuvering  target  -  single  IR  sensor  in  the  presence  of  measurement  noise. 
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Figure  4.15  illustrates  that  the  error  remains  high  during  the  turn,  but  it  is  also 
significantly  reduced  compared  to  all  the  previous  cases.  In  addition,  the  filter  with  the 
two  sensors  recovers  quickly  after  the  turn  and  seeks  the  real  path  of  the  target.  The  chi- 
squared  values  are  much  higher  than  all  the  previous  cases  but  they  last  for  shorter  sam¬ 
ple  periods,  which  means  that  the  tracking  effort  is  large  during  the  turn  and  minimal  at 
all  the  other  times. 

F.  SUMMARY 

We  have  investigated  a  classical  algorithm  applied  to  two  modem  sensors.  It  is 
obvious  that  the  results  are  significantly  improved  in  the  case  where  two  sensors  are  used 
for  tracking.  In  real  applications  more  than  two  sensors  may  be  used  at  the  same  time.  In 
addition,  three  dimensional  cases  are  examined  in  the  literature  [31].  Other  techniques  not 
explored  in  the  present  thesis  may  also  be  used.  Some  alternative  methods  involve  parti¬ 
cle  filters  [29],  and  combined  IR  and  radar  sensors  with  EKF  application  [32]. 
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V.  CONCLUSIONS  AND  FUTURE  WORK. 


We  have  modeled  the  behavior  of  QWIPs  in  various  conditions  of  operation 
through  a  computer  algorithm  and  we  have  simulated  the  behavior  of  an  infrared  sensor 
in  the  bearing-only  tracking  application.  Our  goal,  optimization  of  the  performance  of 
QWIPs  and  their  application  to  target  tracking,  was  achieved  by  implementing  in  a  com¬ 
puter  fundamental  quantum  mechanical  equations  and  nonlinear  control  theory,  respec¬ 
tively.  Our  conclusions  are  summarized  per  chapter  along  with  suggestions  for  future 
work. 

A.  CONCLUSIONS 

In  Chapter  II  we  began  with  the  fundamental  quantum  mechanical  concept,  the 
Schrodinger  equation.  By  application  of  this  equation  in  the  quantum  structures  of  inter¬ 
est  we  studied  their  behavior.  Then,  using  a  well  known  mathematical  technique,  the 
method  of  the  transfer  matrix  (TMM),  we  were  able,  not  only  to  have  our  equations  in  a 
consistent  form,  but  also  to  implement  them  inside  the  simulation  program  and  study 
them.  During  the  program  simulations  we  were  able  to  observe  theoretically  predicted 
and  experimentally  observed  concepts,  such  as  tunneling  through  the  quantum  well  bar¬ 
rier  and  leakage  current.  Finally  we  developed  a  set  of  programming  functions  using  as¬ 
ymptotic  expressions  of  Airy’s  functions  to  study  the  behavior  of  our  structure  in  the  case 
of  very  low  biases. 

Chapter  III  provided  a  discussion  of  the  fabricated  QWIP.  A  single  quantum  well 
of  the  multi-quantum  structure  was  examined  in  various  cases.  Namely,  we  studied  the 
behavior  of  bound,  quasibound  and  continuous  states  in  the  case  of  forward  and  reverse 
bias,  both  above  and  below  the  arbitrarily  placed  threshold  and  in  the  case  of  absence  of 
electric  field.  We  examined  the  concept  of  bound  and  continuous  coefficients  of  the 
wavefunctions  used  in  the  TMM  and  we  developed  their  value  in  every  region  of  the  well 
and  for  the  bias  used.  In  addition,  we  examined  the  perfonnance  of  the  detector  by  study¬ 
ing  its  oscillator  strength  using  characteristic  biases  and  we  found  the  theoretically  pre¬ 
dicted  rapid  oscillations  when  we  apply  very  low  biases.  Moreover,  we  compared  the  ex¬ 
perimental  and  the  simulation  results  and  we  concluded  that  they  match.  Finally,  we  de¬ 
veloped  and  simulated  the  concept  of  photocurrent  in  the  quantum  well. 
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Finally,  Chapter  IV  presented  the  basic  concepts  of  nonlinear  extended  Kalman 
filtering.  We  developed  the  mathematical  model  of  initialization  of  the  filter  using  two 
sets  of  synchronous  bearings.  The  algorithm  was  initialized  by  the  user  and,  based  in  that 
geometry,  various  cases  were  studied,  either  using  one,  or  two  sensors  with  and  without 
the  presence  of  ambient  noise.  Results  and  conclusions  were  made  for  the  optimal  design 
of  a  bearing  only  tracking  system  based  on  IR  sensors. 

B.  FUTURE  WORK 

There  are  many  concepts  based  on  this  work  that  further  research  is  recom¬ 
mended.  For  instance,  using  the  simulation  program  developed,  we  can  design  the  desired 
QWIP,  predicting  its  behavior  in  the  laboratory.  In  addition,  we  can  predict  the  necessary 
bias  to  observe  the  theoretically  predicted  constructive  interference  oscillations.  Further¬ 
more,  we  can  design  the  optimal  quantum  well  to  minimize  the  quantum  tunneling  and  to 
maximize  the  photocurrent  before  testing  it  in  the  laboratory. 

Finally,  it  is  recommended  to  expand  the  EKF  simulation  program,  using  the 
same  initialization  process,  so  that  the  number  of  the  sensors  exceeds  two.  In  that  way  we 
can  simulate  sensors  arrays  and  the  expected  tracking  results  will  be  better.  Furthennore 
we  can  expand  the  simulation  into  three  dimensional  cases.  Concluding,  it  is  strongly 
recommended  to  study  similar  cases  using  particle  filters  and  compare  the  results  with  the 
EKF,  especially  in  three  dimensional  and  in  terrain  mapping  cases  using  synchronous  or 
asynchronous  bearings. 
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APPENDIX  A:  QWIP  SIMULATION  PROGRAMS 


In  the  present  Appendix  we  present  selected  programs  of  the  simulation  of  the 
InGaAs /  AlGaAs  quantum  well.  The  core  of  the  simulation  is  EInGaAs.m.  Once  this 
program  is  executed,  every  other  sub-matlab  file  is  called  automatically. 

%**  Simulation  Of  Performance  Of  Infrared  Photodetectors******* 

%  E_InGaAs 

%LT  Psarakis  Eftychios  Hellenic  Navy 
clc 

clear  all 
global  E  E_c 
%defining  the  constants 
hbar  =  1.05557e-34;  %[J*s] 
mO  =  9.10939e-31;  %[kg] 
q  =  1 .602 1 8e- 19;  %[C] 

F  =  input('Electric  Field  in  [V/m]:  ’); 

%  Percentage  in  region 
x=[0, 0.3, 0.1,0]; 

%  Width  of  the  region 

width=[inf,40,40,inf]; 

ii=4; 

for  tt=  1  :ii 

%X(i)=input('Percentage  of  A1  (0-1)  in  region  '); 

%width(i)  =  input(’Give  width:  ’); 
end 

y  =  0;  %Initial  value  of  y  &  z  to  correct  the  logic  of  the  later  condition 
z  =  0; 
for  tt=2:3 

L(l)  =  width(l)*le-10; 

L(2)  =  width(2)*le-10; 

L(tt+ 1 )  =  (width(tt)+width(tt+ 1  ))*  1  e- 1 0; 
end 

if  F>-2*10A5  &  F<10A5 
%EwellNoBiasAiry 
EwellsmallBiasAiry 
%EwellNoBias 
break 
end 

if  F>=  1 0A5 
for  tt=l:4 
if  width(tt)==inf 
Z_dir(tt)=0; 
else 
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Zdir(tt)  =  L(tt); 
end 
end 
Zdir; 
end 

if  F<=-2*  10A5 
for  tt=l:4 
if  width(tt)==inf 
Z_dir(tt)=0; 
else 

Zdir(l)  =  L(3); 

Z_dir(2)  =  L(2); 
end 
end 

x=[x(4),x(3),x(2),x(l)] 

end 

Zdir 

eg(l)  =  1.519  +  1.247.*x(l); 

eg(2)  =  1.519-  1.102.*x(2); 
eg(3)  =  1.519  -  1.102.*x(3); 
eg(4)  =  1.519  +  1.247  *x(4); 
eg=[eg(l)  eg(2)  eg(3)  eg(4)]; 
for  tt=l:ii 
if  F>0 

V(tt)  =  0.62*(eg(tt)  -  eg(2)); 
else 

V(tt)  =  0.62*(eg(tt)  -  eg(3)); 
end 

m_eff(tt)  =  l/((x(tt)/(0.028*m0))  +  ((l-x(tt))/(0.067*m0))); 

C(tt)  =  ((2*m_eff(tt)*q)/((F*hbar)A2))A(l/3); 
end 

for  tt=l:ii-l 

sigma(tt)  =  ((m_eff(tt)*C(tt+l))/(m_eff(tt+l)*C(tt))); 
end 

if  F  <  0; 

DD  =  0; 

V2_  =  V(l); 
else 

DD  =  -F*L(2); 

V2_  =  V(  1 )-  F*L(3); 
end 

%&&&&&&&&&&&&&&&&&&&&continous  states&&&&&&&&&&&&&&&& 
ac  =1; 
n  =  0; 

11=1; 

%by  varing  dE  we  can  ajust  the  number  of  states  found 


70 


E_c  =min(V(l),V(4))-.03; 

Vb=input('give  Vb  in  eV 
disp('The  continuus  states  energy  are:’) 

%dE_c=  1.2490e-004; 
dE_c=(Vb-E_c)/l  00 1 ; 
while  E_c  <  Vb 
E_c  =  E_c  +  dE_c; 
for  jj=l:4 

eta(jj)  =  (E_c  -  V(jj)); 
end 

for  jj=l:4 

AlphaO'j)  =  -C(jj)*(F*Z_dir(jj)  +  eta(jj)); 

Alphaout(ll,jj)=[Alpha(:  ,jj )] ; 
end 

for  jj=l:3 
alpha(l)=0; 

alpha(jj+l)=  -C(jj + 1 )  *  (F*  Z_dir(jj  )  +  eta(jj+l)); 

alphaout(ll,jj+ 1  )=[alpha(:  ,jj+ 1 )] ; 

end 

forjj=l:3 

ml  l_c=(pi)*(airy(0,alpha(jj+l))*airy(3,Alpha(jj))- 
sigma(jj )  *  (airy(2,  Alpha(jj  ))*  airy  ( 1  ,alpha(jj + 1 )))) ; 

ml2_c=  (pi)*(airy(2,alpha(jj+l))*airy(3,Alpha(jj)) 

sigma(jj )  *  (airy(2,  Alpha(jj  ))*  airy  (3  ,alpha(jj + 1 )))) ; 

m2  l_c=  (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(l,alpha(jj+l))) 

airy(  1  ,Alpha(jj))*airy(0,alpha(jj+l ))); 

m22_c=  (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(3,alpha(jj+l))) 

airy(  1 ,  A I  p ha(j  j ) ) *  ai  ry (2  ,a  1  p h a(j'j  + 1 ))); 

%checking  the  output  data 
ml  lout_c(ll,jj)=[ml  l_c]; 
m  1 2out_c(ll,jj  )=[m  1 2_c] ; 
m2  lout_c(ll,jj)=[m2  l_c]; 
m22out_c(ll,jj  )=[m22_c] ; 

M_c( :  ,jj  )=[m  1 1  _c  ;m  1 2_c;m2 1  _c  ;m22_c] ; 
end 

%finding  Mn 

M12_c=[M_c(l,l),M_c(2,l);M_c(3,l),M_c(4,l)]*[M_c(l,2),M_c(2,2);M_c(3,2), 

M_c(4,2)]; 

M23_c=[M_c(l,2),M_c(2,2);M_c(3,2),M_c(4,2)]*[M_c(l,3),M_c(2,3);M_c(3,3), 

M_c(4,3)]; 

M123_c=[M_c(l,l),M_c(2,l);M_c(3,l),M_c(4,l)]*[M_c(l,2),M_c(2,2);M_c(3,2) 

,M_c(4,2)]*[M_c(1,3),M_c(2,3);M_c(3,3),M_c(4,3)]; 

disp(E_c) 

Eoutcont(ac)  =  E_c; 

Epsicont; 
hold  on 
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r  =  0:.l:width(2)+width(3); 
plot(r,Eout_cont(ac),'r') 
text(0,Eout_cont(ac)+.0 1  ,['E_',num2str(ac+2),’ 
',sprintf('%4.3f,Eout_cont(ac)),'  eV']) 

plot((X*lelO),(real(Psil_cont(ac,:)/5e5)+E_c),'k') 

plot((Y*lelO),(real(Psi2_cont(ac,:)/5e5)+E_c),'k') 

plot((Z*lelO),(real(Psi3_cont(ac,:)/5e5)+E_c),'k') 

plot((T*lelO),(real(Psi4_cont(ac,:)/5e5)+E_c),'k') 

ac=ac+l; 

11=11+1; 

end 

text(-40,1.15*V(l),['V(l)  =  ',sprintf('%4.3f,V(l)),'  eV’]) 
xlabel('X  (Angstroms)') 
ylabel('Potential  Height  (eV)') 
title('Quantum  Well') 

%&&&&&&&&&&&&&&&&&Endofcontinuous  states&&&&&&&&&&& 
% * ********************  *Ewell2** ************************** 

hold  on 

xl  =  -100:. 1:0; 

hi  =  0:.001:V(1); 

plot(x  1 ,( V(  1  )-F*(x  1  *  1  e- 1 0)),'LineWidth',3 ) 
x2  =  0:.l:width(2); 

x3  =  width(2):.l:(width(2)+width(3)); 
if  F  >  0 

h2  =  (-F*F(2)):.001:(V(3)-(F*F(2))); 
h3  =  (V(3)-(FH.(3))):.001:(V(1)-(F*F(3))); 
plot(x3,((V(3)-F*(x3 *  1  e- 1 0)))  ,’LineWidth’,3 ) 
end 

if  F  <  0 

h2  =  (-F*F(2)):.001:(V(2)-F*F(2)); 
h3  =  (V(2)-F*F(3)):.001:(V(1)-F*F(3)); 
plot(x3,((V(2)-F*(x3 *  1  e- 1 0)))  ,’FineWidth’,3 ) 
end 

x4  =  (width(2)+width(3)):.l:(width(2)+width(3)+100); 
plot(x4,((V(  1  )-F*(x4*  1  e- 1 0))),’FineWidth’,3 ) 
plot(x2  ,(-F  *  (x2  *  1  e- 1 0 )) ,  'Fine  W  idth' ,  3 ) 
plot(0,hl  ,'FineWidth’,  1 ) 
plot(width(2),h2,'FineWidth',l) 
plot((width(2)+width(3)),h3,'FineWidth',  1 ) 
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%  %*********************BOTJND  STATES*********************** 

clear  a  n  11  dE  E  eta  alpha  Alpha  E; 
a  =1; 
n  =  0; 

u=i; 

dE  =  5*10A-4; 

E  =  DD; 

%Evaluate  the  energy  from  the  bottom  to  the  top  of  the  well; 
while  E  <  V2_-dE 
if  E  >  V2_-10A-3 
dE  =  10A-6; 
end 

E  =  E  +  dE; 
for  tt=  1 :4 

eta(tt)  =  (E  -  V(tt)); 
end 

for  jj=l:4 

AlphaO'j)  =  -  C  0  J )  *  (F  *  Zdir  O'j )  +  eta(jj)); 
end 

for  jj=l:3 

alpha(jj+l)=  -C(jj+ 1 )  *  ( F  *  Zdi  r(jj )  +  eta(jj+l)); 
end 

for  jj=l:3 

ml  1=  (pi)*(airy(0,alpha(jj+l))*airy(3,Alpha(jj)) 

sigma  (jj  )*(airy(2,Alpha(jj))*airy(  1  ,alpha(jj+l )))); 

m  1 2=  (p  i)  *  (airy  (2  ,alpha(jj  + 1 ) )  *  airy  (3 ,  Alpha  (jj ) ) 

si gma(jj )  *  (airy (2 ,  A lph a(jj ))  *  ai  ry (3 , alph a(jj + 1 )))  ) ; 

m2 1=  (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(  1  ,alpha(jj+l ))) 

airy(  1 ,  AlphaO'j  ))*  airy  (0  ,alpha(jj + 1 ))); 

m22=  (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(3,alpha(jj+l))) 

airy(  1 ,  Alpha(jj  ))*airy(2  ,alpha(jj+ 1 ))); 

%checking  the  output  data 
ml  lout(ll,jj)=[ml  1]; 
m  1 2out(ll,jj)=[m  12]; 
m2 1  out(ll,jj)=[m2 1  ] ; 
m22out(ll,jj)=[m22]; 

M(:,jj)=[ml  I;ml2;m21;m22]; 
end 

%finding  Mn 

M 1 2=[M(  1 , 1  ),M(2, 1  );M(3 , 1  ),M(4, 1 )]  *  [M(  1 ,2),M(2,2);M(3 ,2),M(4,2)] ; 
M23=[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,3),M(2,3);M(3,3),M(4,3)]; 

M  123=[M(  1 , 1  ),M(2, 1  );M(3, 1  ),M(4,1  )]*  [M(  1 ,2),M(2,2);M(3,2),M(4,2)]*[M(  1 ,3), 
M(2,3);M(3,3),M(4,3)]; 
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M22  =  real(M123(2,2));  %in  case  we  have  three  matrices 
%checking  the  output  data 
M22out(ll)=[M22] ; 

EM22(11)=[E]; 

%%%%%%%end  of  check%%%%%%%%%%%%%%% 


z  =  y; 
y  =  M22; 

11=11+1; 

if  z*y  <  0  %Detennine  the  approximate  eigen  energy  first 
n  =  n  +1; 

if  n  <  2  %Increase  the  acuaracy  of  eigen  energy 
E  =  E-dE; 
dE  =  le-6; 
y  =  z; 
end 
end 
if  n  >  1 

disp('The  ground  state  energy  is:’) 
disp(E) 

Eout(a)  =  E;  %Stores  energy  levels  into  an  array 
Epsibound; 
dE  =  5e-4; 
n  =  0; 
hold  on 

r  =  0:.l:width(2)+(a-l)*width(3); 
if  a  <  3 

plot(r,Eout(a),’r’) 

text(0,Eout(a)+.01,['E_',num2str(a),'  =  ',sprintf('%4.3f,Eout(a)),'  eV']) 
plot((X*  1  e  1 0),(real(Psi  l_bound(a,:)/5e5  )+Eout(a)),'k') 
plot((Y*lelO),(real(Psi2_bound(a,:)/5e5)+Eout(a)),'k') 
plot((Z*lelO),(real(Psi3_bound(a,:)/5e5)+Eout(a)),'k') 
plot((T*lelO),(real(Psi4_bound(a,:)/5e5)+Eout(a)),'k') 
end 

a  =  a+ 1 ; 
end 
end 
ZETA; 

o/o strength* 
for  ai=l:ac-l; 

DE(ai)=(Eout_cont(ai)-Eout); 

f_oscil(ai)=2*(m_eff(l))/(hbarA2).*DE(ai)*q.*Zeta_final(ai)'; 

end 

figure 

plot(DE,real(f_oscil)) 
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xlabel('YDeltaE') 
ylabel('oscillator  strength  f ) 
title('oscillator  strength  Vs  YDeltaE') 
for  ai=l:ac-l; 

Prod_f_oscil_DOS(ai)=f_oscil(ai).*DOS(ai)*q; 

OSCSTR(ai)=f_oscil(ai).*DOS(ai)*q*dE_c; 

OSCSTR_INTEG=sum(OSCSTR); 

end 

figure 

plot(DE,real(Prod_f_oscil_DOS)) 

xlabel('\DeltaE') 

ylabel('(oscillator  strength)*(density  of  states)') 
title('oscillator  strength  multiplied  by  density  of  states  Vs  YDeltaE') 
disp('The  (oscillator  strength)*(density  of  states)  integral  is') 
round(OSCSTRINTEG) 


%  *******************  * *  Absorption* ********************** 


q=  1.602 18e- 19;  %[C] 
h=6.626*10A-34; 
m_eff(l)=1.0e-03 1*0. 6103; 

N_d=le24; 

e0=8.85e-12; 

n_r=3.5; 

c=3*10A8; 

Cb_c=pi*(N_d*q*hbarA2)./(2*m_eff(  1  ).A2*eO*n_r*c); 
total_length=3  80*  1 0A- 10; 
for  kk=l:ac-l 


beta4=((hbarA2)/(2*m_eff(4)*F))A(l/3); 

%region  IV 

quantum_eff=total_length*Cb_c .  *Prod_f_oscil_DOS ; 

J4R(kk)=hbar/(4*m_eff(4)*pi*abs(beta4))*(b3_cont(kk).*A4(kk)+b3_cont(kk).*( 

i*B4(kk))).*conj((b3_cont(kk).*A4(kk)+b3_cont(kk).*(i*B4(kk))))*q; 

%  I=J_total*Area; 
end 

J_total=quantum_eff.  *  J4R' .  *  q; 
figure 

plot((h*c)./(q*DE),real(J_total)/max(real(J_total)),'b’) 
hold  on 

plot((h*c)./(q*DE),real(Cb_c*Prod_f_oscil_DOS)./max(real(Cb_c*Prod_f_oscil_ 

DOS)),’r’) 

xlabel('Wavelength-\mum') 

ylabel('Nonnalized  photocurrent,  Nonnalized  absorption') 
legendCphotocurrent', 'absorption') 
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%**Simulation  Of  Performance  Of  Infrared  Photodetectors'1”1”1”1”15** 
%  EwellNoBias 

%LT  Psarakis  Eftychios  Hellenic  Navy 
hbar  =  1.05557e-34;  %[J*s] 
mO  =  9.10939e-31;  %[kg] 
q=  1.602 18e- 19;  %[C] 

F  =F; 

%  Percentage  in  region 
x=[0, 0.3, 0.1,0]; 

%  Width  of  the  region 

width=[inf,40,40,inf]; 

ii=4; 

Y  =  0; 

z  =  0; 

for  tt=2:ii-l 

L(l)  =  width(l)*le-10; 

L(2)  =  width(2)*le-10; 

L(tt+1)  =  (width(tt)+width(tt+l))*le-10; 
end 

eg(l)  =  1.519  +  1.247.*x(l); 
eg(2)=  1.519-  1.102.*x(2); 
eg(3)  =  1.519-  1.102.*x(3); 
eg(4)=  1.519+  1.247.*x(4); 
eg=[eg(l)  eg(2)  eg(3)  eg(4)]; 

V(l)=0.62*(eg(l)-eg(2)); 

V(2)=0; 

V(3)=0.62*(eg(3)-eg(2)); 

V  (4)=0.62*(eg(4)-eg(2)); 

V=[V(1),V(2),Y(3),V(4)]; 
for  jj=l:4 

m_eff(jj)  =  1  /  ((x(jj  )/(0 . 02 8  *  m0))  +  ((l-x(jj))/(0.067*m0))); 
end 

for  t=l:ii 

if  width(t)==inf 
Z_dir(t)=0; 
else 

Z_dir(t)  =  L(t); 
end 
end 
ac  =1; 
n  =  0; 

n=i; 

%by  varing  dE  we  can  ajust  the  number  of  states  found 
E  =max(  V  ( 1 ),  V  (4)); 


76 


jj=0; 

Vb=input('give  Vb  in  eV  :') 
dE_c=10A-4; 

disp('The  continuus  states  energy  are:’) 
while  E<Vb 
E=E+dE_c; 
for  nn=l:4 

k(nn)=sqrt((2  *  m_eff(nn)  *(E-V(nn))*  q)/hbar  A2 ) ; 
end 

for  nn=l  :3 

XI(nn)=(k(nn)  *m_eff(nn+ 1 ))/ (k(nn+ 1 )  *m_eff(nn)) ; 
end 

for  jj=l:3 

ml  l_c=  1  /2*((  1  +XI(jj  ))*exp(i*(k(jj  )-k(jj+ 1 ))  *Z_dir(j  j ))) ; 
ml2_c=  1  /2  *  ( ( 1  -XI  (j j  )  )  *  exp  ( -  i  *  (k  (jj ) + k  (jj  + 1 ) )  *  Z_d  i  r(jj ) ) ) ; 
m2 1  _c=  1  /2  *  (( 1  -XI(jj ))  *  exp(i  *  (k(jj  )+k(jj + 1 ) )  *Z_dir(jj  ))) ; 
m22_c=  l/2*((l+XI(jj))*exp(-i*(k(jj)-k(jj+l))*Z_dir(jj))); 

M_c( :  ,jj  )=[m  1 1  _c  ;m  1 2_c  ;m2 1  _c  ;m22_c] ; 
ml  lout_c(ll,jj)=[ml  l_c]; 
m  1 2  out_c(ll,jj  )=[m  1 2_c] ; 
m2 1  out_c(ll,jj)=[m2  l_c] ; 
m22out_c(ll,jj)=[m22_c] ; 
end 

%finding  Mn 

M2  l_c=[M_c(  1 ,2),M_c(2,2);M_c(3,2),M_c(4,2)]  *  [M_c(  1 , 1  ),M_c(2, 1  );M_c(3, 1 ), 
M_c(4,l)]; 

M32_c=[M_c(1,3),M_c(2,3);M_c(3,3),M_c(4,3)]*[M_c(1,2),M_c(2,2);M_c(3,2), 

M_c(4,2)]; 

M321_c=[M_c(l,3),M_c(2,3);M_c(3,3),M_c(4,3)]*[M_c(l,2),M_c(2,2);M_c(3,2) 
,M_c(4,2)]  *  [M_c(  1 , 1  ),M_c(2, 1  );M_c(3 , 1  ),M_c(4, 1 )] ; 

M22_c  =real(M32 1 _ c(2,2)); 

%checking  the  output  data 
disp(E) 

Eoutcont(ac)  =  E; 
integrationcont; 
hold  on 

r  =  0:.l:width(2)+width(3); 
plot(r,Eout_cont(ac),'r’) 

text(0,Eou_contt(ac)+.0 1  ,['E_',num2str(ac),'  = 

',sprintf('%4.3f,Eout_cont(ac)),'  eV']) 

plot(reg  1  *  1 0A 1 0,(real(psi  l_cont(ac)/5e5  )+Eout_cont(ac)),'k’) 
plot(reg2*10A10,(real(psi2_cont(ac)/5e5)+Eout_cont(ac)),'k') 
plot(reg3  *  1 0A 1 0,(real(psi3_cont(ac)/5e5)+Eout_cont(ac)),'k') 

%  plot(reg4*10A10,(real(psi4_cont(ac)/5e5)+Eout_cont(ac)),'k') 
ac=ac+l; 
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end 

xlabel('X  (Angstroms)') 
ylabel('Potential  Height  (eV)') 
title('Quantum  Well') 

/o  IV  ell 

hold  on 

xl  =  -100:. 1:0; 

hi  =  0:.001:V(1); 

plot(xl ,(V(  1  )-F*(xl*  1  e- 1 0)))  %,'LineWidth',5) 
x2  =  0:.l:width(2); 

x3  =  width(2):.l:(width(2)+width(3)); 
h2  =  (-F*L(2)):.001:(V(3)-(F*L(2))); 
h3  =  (V(3)-(F*L(3))):.001:(V(1)-(F*L(3))); 
plot(x3 ,((V(3)-F*(x3 *  1  e- 1 0))))  %,’LineWidth’,5) 
x4  =  (width(2)+width(3)):.l:(width(2)+width(3)+100); 
plot(x4,((V(  1  )-F*(x4*  1  e- 1 0))))  %,’LineWidth’,5) 
x5  =  0:.l:(width(2)+width(3)); 
p  lot(x2  ,(-F  *  (x2  *  1  e- 1 0 ) ))  % ,  'Line Width' ,  5  ) 
plot(0,hl)  %,'LineWidth',8) 
plot(width(2),h2)  %,’LineWidth',3) 
plot((width(2)+width(3)),h3)  %,'LineWidth',5) 

a  =1; 
n  =  0; 

11=1; 

dE  =  5*10A-4; 

E  =0; 
jj=0; 

while  E<max(V) 
if  E  >  max(V)-10A-3 
dE  =  10A-6; 
end 

E=E+dE; 
for  nn=l:4 

k(nn)=sqrt((2  *m_eff(nn)  *(E- V(nn))  *  q)/hbarA2) ; 
end 

for  nn=l:3 

XI(nn)=(k(nn)*m_eff(nn+ 1  ))/(k(nn+ 1  )*m_eff(nn)); 
end 

for  jj=l:3 

ml  1=  l/2*((l+XI(jj))*exp(i*(k(jj)-k(jj+l))*Z_dir(jj))); 
ml  2=  l/2*((l-XI(jj))*exp(-i*(k(jj)+k(jj+l))*Z_dir(jj))); 
m21=  l/2*((l-XI(jj))*exp(i*(k(jj)+k(jj+l))*Z_dir(jj))); 


78 


m22=  l/2*((l+XI(jj))*exp(-i*(k(jj)-k(jj+l))*Z_dir(jj))); 

M(:,jj)=[ml  I;ml2;m21;m22]; 

ml  lout(ll,jj)=[ml  1]; 

m  1 2out(ll,jj)=[m  12]; 

m2 1  out(ll,jj)=[m2 1  ]; 

m22out(ll,jj)=[m22]; 

end 

%finding  Mn 

M2 1=[M(  1 ,2),M(2,2);M(3,2),M(4,2)]*[M(  1 , 1  ),M(2, 1);M(3, 1  ),M(4, 1 )]; 
M32=[M(1,3),M(2,3);M(3,3),M(4,3)]*[M(1,2),M(2,2);M(3,2),M(4,2)]; 

M32 1 =[M(  1 ,3),M(2,3);M(3,3),M(4,3)]*  [M(  1 ,2),M(2,2);M(3,2),M(4,2)]  *  [M(  1,1), 
M(2,1);M(3,1),M(4,1)]; 

M22  =real(M32 1(2,2)); 

%checking  the  output  data 
M22out(ll)=[M22]; 

EM22(11)=[E]; 

%%%%%%%end  of  check%%%%%%%%%%%%%%% 

z=y; 

y=M22  ; 

11=11+1; 
if  z.*y  <  0 
n  =  n  +1; 

ifn  <  2  %Increase  the  acuaracy  of  eigen  energy 
E  =  E-dE; 
dE  =  le-6; 
y  =  z; 
end 
end 
if  n  >  1 

disp(’The  ground  state  energy  is:') 
disp(E) 

Eout(a)  =  E; 
xx=E; 
integration; 
dE  =  5e-4; 
n  =  0; 
hold  on 

r  =  0:.l:width(2)+(a-l)*width(3); 

Eout(a)  =  E; 
hold  on 
if  a<3 

plot(r,Eout(a),'r') 

text(0,Eout(a)+.01,['E_',num2str(a),'  =  ',sprintf('%4.3f,Eout(a)),'  eV']) 
plot(reg  1  *  1  e  1 0,(real(Psi  l_bound/5e5)+E),'k') 
plot(reg2  *  1  e  1 0,(real(Psi2_bound/5e5)+E),'k') 
plot(reg3  *  1  e  1 0,(real(Psi3_bound/5e5)+E),'k') 
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plot(reg4*  1  e  1 0,(real(Psi4_bound/5e5)+E),'k') 
end 
a=a+ 1 ; 
end 
end 

ZETA1; 

%******************************************************** 
% **** **********osciiiator  strength*********************** 

for  ai=l:ac-l; 

DE(ai)=(Eout_cont(ai)-Eout(l)); 

f_oscil(ai)=2*(m_eff(l))/(hbarA2).*DE(ai)*q.*Zeta_final(ai)'; 

end 

figure(2) 

plot(DE,real(f_oscil)) 

xlabel('\DeltaE') 
ylabel('oscillator  strength  f ) 
title('oscillator  strength  Vs  YDeltaE') 
for  ai=l:ac-l; 

Prod_f_oscil_DOS(ai)=f_oscil(ai).*DOS(ai)*q; 

OSCSTR(ai)=f_oscil(ai).*DOS(ai)*q*dE_c; 

OSCSTR_INTEG=sum(OSCSTR); 

end 

figure(3) 

plot(DE,real(Prod_f_oscil_DOS)) 

xlabel('\DeltaE') 

ylabel('(oscillator  strcngth)*(dcnsity  of  states)') 
title('oscillator  strength  multiplied  by  density  of  states  Vs  YDeltaE') 
disp('The  (oscillator  strength)* (density  of  states)  integral  is') 
round(OSCSTRINTEG) 


%**Simulation  Of  Performance  Of  Infrared  Photodetectors******* 
%  integration 

%LT  Psarakis  Eftychios  Hellenic  Navy 
dx  =  10A- 1 1 ; 

regl=  -300e-10:dx:0; 

reg2  =  0:dx:L(l,2); 

reg3  =  L(l,2):dx:L(l,3); 

reg4  =  L(l,3):dx:(L(l,3)+300e-10); 


A1=0; 

Bl=l; 

A2=M(2,1); 
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B2=M(4,1); 

A3=M2 1(1,2); 

B3=M2 1(2,2); 

A4=M32 1(1,2); 

B4=0; 

psil_bound=Al*exp(i*k(l)*regl)+Bl*exp(-i*k(l)*regl); 

psi2_bound=A2*exp(i*k(2)*reg2)+B2*exp(-i*k(2)*reg2); 

psi3_bound=A3*exp(i*k(3)*reg3)+B3*exp(-i*k(3)*reg3); 

psi4_bound=A4*exp(i*k(4)*reg4)+B4*exp(-i*k(4)*reg4); 

integrall=sum((psil_bound.*conj(psil_bound))*dx); 

integral2=sum((psi2_bound.*conj(psi2_bound))*dx); 

integral3=sum((psi3_bound.*conj(psi3_bound))*dx); 

integral4=sum((psi4_bound.*conj(psi4_bound))*dx); 

total_integral=integral  1  +integral2+integral3+integral4; 

b  1  =(sqrt(total_integral)); 

Psi  1  _bound=  1/bl  *psil  bound; 

Psi2_bound=l/b  1  *psi2_bound; 

Psi3_bound=l/bl*psi3_bound; 

Psi4_bound=l/b  1  *psi4_bound; 

Integral  1  =sum((Psi  1  bound.  *  conj  (Psi  1  bound))  *  dx); 
Integral2=sum((Psi2_bound.*conj(Psi2_bound))*dx); 
Integral3=sum((Psi3_bound.*conj(Psi3_bound))*dx); 
Integral4=sum((Psi4_bound.*conj(Psi4_bound))*dx); 

T  otal_integral=Integral  1  +Integral2+Integral3+Integral4; 

%**Simulation  Of  Performance  Of  Infrared  Photodetectors******* 

%  integration_cont 

%LT  Psarakis  Eftychios  Hellenic  Navy 
dx  =  10A-1 1; 

regl=  -300e-10:dx:0; 

reg2  =  0:dx:L(l,2); 

reg3  =  L(l,2):dx:L(l,3); 

reg4  =  L(l,3):dx:(L(l,3)+300e-10); 

Al_c(ac)=l/sqrt(2*pi); 

Bl_c(ac)=-Al_c(ac)*M321_c(2,l)/M321_c(2,2); 

A2_c(ac)=M_c(  1 ,  1)*  A  l_c(ac)+M_c(2, 1)*B  l_c(ac); 
B2_c(ac)=M_c(3,l)*Al_c(ac)+M_c(4,l)*Bl_c(ac); 

A3_c(ac)=M2  l_c(  1 , 1)*  A  l_c(ac)+M2  l_c(  1 ,2)*B  l_c(ac); 
B3_c(ac)=M21_c(2,l)*Al_c(ac)+M21_c(2,2)*Bl_c(ac); 

A4_c(ac)=M32  l_c(  1  ,l)*Al_c(ac)+M32  l_c(  1 ,2)*B  l_c(ac); 

B4_c(ac)=0; 

psil_cont(ac,:)=Al_c(ac)*exp(i*k(l)*regl)+Bl_c(ac)*exp(-i*k(l)*regl); 

psi2_cont(ac,:)=A2_c(ac)*exp(i*k(2)*reg2)+B2_c(ac)*exp(-i*k(2)*reg2); 

psi3_cont(ac,:)=A3_c(ac)*exp(i*k(3)*reg3)+B3_c(ac)*exp(-i*k(3)*reg3); 

psi4_cont(ac,:)=A4_c(ac)*exp(i*k(4)*reg4)+B4_c(ac)*exp(-i*k(4)*reg4); 
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%**Simulation  Of  Performance  Of  Infrared  Photodetectors******* 

%  ZETA1 

%LT  Psarakis  Eftychios  Hellenic  Navy 
dx  =  10A-1 1; 

regl=  -300e-10:dx:0; 
reg2  =  0:dx:L(l,2); 
reg3  =  L(l,2):dx:L(l,3); 
reg4  =  L(l,3):dx:(L(l,3)+300e-10); 
dim=max(reg4)-min(reg  1 ); 
for  kk=l:ac-l 

Zital(kk,:)=conj(Psil_bound(l,:)).*regl.*psil_cont(kk,:).*dx; 
Zita2(kk,:)=conj‘(Psi2_bound(l,:)).*reg2.*psi2_cont(kk,:).*dx; 
Zita3(kk,:)=conj(Psi3_bound(l,:)).*reg3.*psi3_cont(kk,:).*dx; 
Zita4(kk,:)=conj(Psi4_bound(l,:)).*reg4.*psi4_cont(kk,:).*dx; 
OSC_integrall(kk,:)  =  sum(Zital(kk,:)); 

OSC_integral2(kk,:)  =  sum(Zita2(kk,:)); 

OSC_integral3(kk,:)  =  sum(Zita3(kk,:)); 

OSC_integral4(kk,:)  =  sum(Zita4(kk,:)); 

O  S  C_flnali(kk, :  )=(0  S  Cintegral  1  (kk, : )+  O  SC_integral2  (kk, : )+ 

OSC_integral3(kk,:)+OSC_integral4(kk,:)); 

Zeta_final(kk,:)=(OSC_finab(kk,:)).*conj(OSC_finali(kk,:)); 

%DOS(kk)=dim./(hbar*pi)*sqrt(m_eff(l)./(2*q*(Eout_cont(kk)- 

Eout_cont(  1  )+0 .00 1 ))); 

DOS(kk)=l/(hbar*pi)*sqrt(m_eff(l)./(2*q*(Eout_cont(kk)-V(l)+0.0001))); 

end 


%**Simulation  Of  Performance  Of  Infrared  Photodetectors******* 
%E  wellsmallB  ias  Airy 
%LT  Psarakis  Eftychios  Hellenic  Navy 
if  F>0 
for  tt=l:4 
if  width(tt)==inf 
Z_dir(tt)=0; 
else 

Zdir(tt)  =  L(tt); 
end 
end 
Zdir 
else 

for  tt=l:4 
if  width(tt)==inf 
Z_dir(tt)=0; 
else 
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Z  dir(l)  =  L(3); 

Z_dir(2)  =  L(2); 
end 
end 

x=[x(4),x(3  ),x(2),x(  1 )] ; 
end 

eg(l)  =  1.519  +  1.247.*x(l); 

eg(2)  =  1.519-  1.102.  *x(2); 
eg(3)=  1.519-  1.102. *x(3); 
eg(4)=  1.519+  1.247.*x(4); 
eg=[eg(l)  eg(2)  eg(3)  eg(4)]; 
for  tt=l:ii 
if  F>0 

V(tt)  =  0.62*(eg(tt)  -  eg(2)); 
else 

V(tt)  =  0.62*(eg(tt)  -  eg(3)); 
end 

m_eff(tt)  =  l/((x(tt)/(0.028*m0))  +  ( (1  -  x  ( U ) )/( 0 . 0 6 7 *  m  ()))); 

C(tt)  =  (((2  *m_eff(tt)  *  q)/((F  *hbar)A2))A(  1/3)); 

end 

for  tt=l:ii-l 

sigma(tt)  =  ((m_eff(tt)*C(tt+l))/(m_eff(tt+l)*C(tt))); 
end 

if  F  <  0; 

DD  =  0; 

V2_  =  V(l); 
else 

DD  =  -F*L(2); 

V2_  =  V(l)-  F*L(3); 
end 
ac  =1; 
n  =  0; 

11=1; 

%by  varing  dE  we  can  ajust  the  number  of  states  found 
E_c  =min(V(l),V(4))-.03; 

Vb=input('give  Vb  in  eV 

dispCThe  continuus  states  energy  are:') 

dE_c  =  (Vb-E_c)/1001; 

while  E_c  <  Vb 

E_c  =  E_c  +  dE_c; 

for  tt=l:4 

eta(tt)  =  (E_c  -  V(tt)); 
end 

for  jj=l:4 

Alpha(jj)  =  -C(jj)*(F*Z_dir(jj)  +  eta(jj)); 
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Alphaout(ll,jj)=[Alpha(:  ,jj )] ; 
end 

for  jj=l:3 

alpha(jj+l)=  -C(jj + 1 )  *  (F*  Z_dir(jj )  +  eta(jj+l)); 

alphaout(ll,jj  )=[alpha( :  ,jj )] ; 

end 

forjj=l:3 

mll_c=  (pi)*(airy(0,alpha(jj+l))*airy(3,Alpha(jj)) 

sigma(jj)*(airy(2,Alpha(jj))*airy(  1  ,alpha(jj+ 1 )))); 

ml2_c=  (pi)*(airy(2,alpha(jj+l))*airy(3,Alpha(jj)) 

sigma(jj )  *  (airy(2 ,  Alpha(jj  ))*  airy  (3  ,alpha(jj + 1 )))) ; 

m21_c=  (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(l,alpha(jj+l))) 

airy(  1  ,Alpha(jj))*airy(0,alpha(jj+l))); 

m22_c=  (pi)*(sigma ( jj )  *  ( a  i  ry  ( 0 ,  A 1  p  h  a ( jj )  )  *  airy  ( 3 ,  alphaj j + 1 )  )) 

airy(  1 ,  Alpi  ha(j'j ))*  airy  (2  ,alph  a(jj  + 1 ))); 

%checking  the  output  data 
ml  lout_c(ll,jj)=[ml  l_c]; 
m  1 2out_c(ll,jj  )=[m  1 2_c] ; 
m2 1  out_c(ll,jj)=[m2  l_c] ; 
m22out_c(ll,jj  )=[m22_c] ; 

M_c( :  jj  )=[m  1 1  _c  ;m  1 2_c;m2 1  _c  ;m22_c] ; 
end 

%  finding  Mn 

M12_c=[M_c(  1 , 1  ),M_c(2, 1  );M_c(3 , 1  ),M_c(4, 1 )]  *  [M_c(  1 ,2),M_c(2,2);M_c(3,2), 
M_c(4,2)]; 

M23_c=[M_c(1,2),M_c(2,2);M_c(3,2),M_c(4,2)]*[M_c(1,3),M_c(2,3);M_c(3,3), 

M_c(4,3)]; 

M123_c=[M_c(l,l),M_c(2,l);M_c(3,l),M_c(4,l)]*[M_c(l,2),M_c(2,2);M_c(3,2) 

,M_c(4,2)]*[M_c(l,3),M_c(2,3);M_c(3,3),M_c(4,3)]; 

disp(Ec) 

Eoutcont(ac)  =  E_c; 

Epsicont; 
hold  on 

r  =  0:.l:width(2)+width(3); 
plot(r,Eout_cont(ac),’rj 

text(0,Eout_cont(ac)+.0 1  ,[’E_',num2str(ac+2),’  = 

’,sprintf('%4.3f,Eout_cont(ac)),'  eV']) 

plot((X*lelO),(real(psil_cont(ac,:)/5e5)+Eout_cont(ac)),’kj 

plot((Y*lelO),(real(psi2_cont(ac,:)/5e5)+Eout_cont(ac)),’kj 

plot((Z*lel0),(real(psi3_cont(ac,:)/5e5)+Eout_cont(ac)),’kj 

plot((T*lelO),(real(psi4_cont(ac,:)/5e5)+Eout_cont(ac)),’kj 

ac=ac+l; 

11=11+1; 

end 

text(-40,1.15*V(l),[’V(l)  =  ',sprintf(’%4.3f , V(  1 )),'  eV’]) 
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xlabel('X  (Angstroms)') 
ylabel('Potential  Height  (eV)') 
title('Quantum  Well') 

hold  on 

xl  =  -100:.  1:0; 

hi  =  0:.001:V(1); 

plot(x  1 ,( V(  1  )-F*(x  1  *  1  e- 1 0)),'LineWidth',3) 
x2  =  0:.l:width(2); 

x3  =  width(2):.l:(width(2)+width(3)); 
if  F  >  0 

h2  =  (-F*L(2)):.001:(V(3)-(F*L(2))); 
h3  =  (V(3)-(F*L(3))):.001:(V(1)-(F*L(3))); 
plot(x3,((V(3)-F*(x3*  le-10)))  ,’LineWidth’,3) 
end 

if  F  <  0 

h2  =  (-F*L(2)):.001:(V(2)-F*L(2)); 
h3  =  (V(2)-F*L(3)):.001:(V(1)-F*L(3)); 
plot(x3,((V(2)-F*(x3*  le-10)))  ,’LineWidth’,3) 
end 

x4  =  (width(2)+width(3)):.l:(width(2)+width(3)+100); 
plot(x4,((V(l)-F*(x4*le-10))),'LineWidth',3) 
plot(x2  ,(-F  *  (x2  *  1  e- 1 0 )) ,  'Line  W  idth' ,  3 ) 
p  lot(0  ,h  1 ,  ’Line  W  idth’ ,  1 ) 
plot(width(2),h2,’LineWidth',  1 ) 
plot((width(2)+width(3)),h3,'LineWidth’,  1 ) 

o^******************  *QQpj]\fp)  STATES  ****************************** 

%  clear  M22  M22out  n  11  dE  E  ml  1  ml 2  m21  m22  alpha  alphaout 
%  clear  ml  lout  ml2out  m2  lout  m22out  M12  M23  M123  Alpha  Alphaout 
ab  =1; 
n  =  0; 

11=1; 

dE  =  5*10A-4; 

E  =  DD; 

%Evaluate  the  energy  from  the  bottom  to  the  top  of  the  well; 
while  E  <  V2_-dE 

if  E  >  V2_-10A-3 
dE  =  10A-6; 
end 

E  =  E  +  dE; 
for  tt=l:4 

eta(tt)  =  (E  -  V(tt)); 
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end 

for  jj=l:4 

Alpha(jj)  =  -C(jj)*(F*Z_dir(jj)  +  eta(jj)); 

Alphaout(ll,jj)=[Alpha(:  ,jj)]; 
end 

for  jj=l:3 

alpha(jj+l)=  -C(jj + 1 )  *  (F*  Z_dir(jj )  +  eta(jj+l)); 

alphaout(ll,jj  )=[alpha( :  ,jj )] ; 

end 

forjj=l:3 

m  1 1 =(pi)  *  (asymptotic  Airy(0,alpha(jj+ 1  ))*  asymptotic  Airy(3  ,Alpha(jj )) 
sigma (jj )  *(asymptotic  Airy(2 ,  Alpha(jj ))  *  asymptotic  Airy(  1  ,alpha(jj + 1 )))); 

ml2=(pi)*(asymptoticAiry(2,alpha(jj+l))*asymptoticAiry(3,Alpha(jj)) 

sigma(jj)*(asymptoticAiry(2,Alpha(jj))*asymptoticAiry(3,alpha(jj+l)))); 

m21= 

(pi)*(sigma(jj)*(asymptoticAiry(0,Alpha(]j))*asymptoticAiry(l,alpha(jj+l)))  -  asymp 
toticAiry(  1  ,Alpha(jj))*asymptoticAiry(0,alpha(jj+ 1 ))); 
m22= 

(pi)  *  (sigma  (jj )  *  (asymptotic  Airy  ( 0 ,  Alpha(jj ) )  *  asymptotic  Airy  (3 ,  alpha(jj + 1 )))  -  asymp 
toticAiry(  1  ,Alpha(jj))*asymptoticAiry(2,alpha(jj+ 1))); 

%checking  the  output  data 
ml  lout(ll,jj)=[ml  1]; 
m  1 2out(ll,jj)=[m  12]; 
m2 1  out(ll,jj)=[m2 1  ] ; 
m22out(ll,jj)=[m22]; 

M(:  ,jj)=  1 0A-40*  [ml  1  ;m  1 2;m2 1  ;m22] ; 
end 

%finding  Mn 

M12=[M(1,1),M(2,1);M(3,1),M(4,1)]*[M(1,2),M(2,2);M(3,2),M(4,2)];  %10A60 
M23=[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,3),M(2,3);M(3,3),M(4,3)];  %10a60 
M 1 23=[M(  1 , 1  ),M(2, 1  );M(3 , 1  ),M(4, 1  )]*  [M(  1 ,2),M(2,2);M(3,2),M(4,2)]  *  [M(  1 ,3), 
M(2,3);M(3,3),M(4,3)];  %10a90 
M22  =  real(M123(2,2)); 

%checking  the  output  data 
M22out(ll)=[M22] ; 

EM22(11)=[E]; 

%%%%%%%end  of  check%%%%%%%%%%%%%%% 

11=11+1; 

z  =  y; 
y  =  M22; 

if  z*y  <  0  %Detennine  the  approximate  eigen  energy  first 

n  =  n  +1; 

ifn  <  2  %Increase  the  acuaracy  of  eigen  energy 
E  =  E-dE; 
dE  =  le-6; 
y  =  z; 
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end 
end 
if  n  >  1 

disp(’The  ground  state  energy  is:') 
disp(E) 

%  disp(M123) 

%  disp(M23) 

Eout(ab)  =  E;  %Stores  energy  levels  into  an  array 
Epsiboundl; 
dE  =  5e-4; 
n  =  0; 
hold  on 

r  =  0:.l:width(2)+(ab-l)*width(3); 
if  ab  <  3 

plot(r,Eout(ab),’r’) 

text(0,Eout(ab)+.01,[’E_',num2str(ab),'  =  ’,sprintf('%4.3f,Eout(ab)),'  eV’]) 

plot((X*lelO),(real(Psil_bound(ab,:)/5e5)+E),'k') 

plot((Y*lelO),(real(Psi2_bound(ab,:)/5e5)+E),'k') 

plot((Z*lel0),(real(Psi3_bound(ab,:)/5e5)+E),'k') 

plot((T*lelO),(real(Psi4_bound(ab,:)/5e5)+E),'k') 

end 


ab  =  ab+1; 
end 
end 
ZETA; 

for  ai=l:ac-l; 

DE(ai)=(Eout_cont(ai)-Eout); 

f_oscil(ai)=2*(m_eff(l))/(hbarA2).*DE(ai)*q.*Zeta_final(ai)'; 

end 

figure(2) 

plot(DE,real(f_oscil)) 


xlabel('\DeltaE') 
ylabel('oscillator  strength  f ) 
title('oscillator  strength  Vs  \DeltaE') 
for  ai=l:ac-l; 

Prod_f_oscil_DOS(ai)=f_oscil(ai).*DOS(ai)*q; 

OSCSTR(ai)=f_oscil(ai).*DOS(ai)*q*dE_c; 

OSCSTR_INTEG=sum(OSCSTR); 

end 

figure(3) 

plot(DE,real(Prod_f_oscil_DOS)) 
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xlabel('\DeltaE') 

ylabel('(oscillator  strength)* (density  of  states)') 

title('oscillator  strength  multiplied  by  density  of  states  Vs  YDeltaE') 

disp(’The  (oscillator  strength)* (density  of  states)  integral  is') 
round(OSCSTRINTEG) 


%**  Simulation  Of  Performance  Of  Infrared  Photodetectors******* 
%Epsi_boundl 

%LT  Psarakis  Eftychios  Hellenic  Navy 
warning  off 
deltax  =  1 0A- 1 1; 

B4=  1; 

A4  =  0; 

A3(ab)  =  10A40*M(2,3)*B4; 

B3(ab)  =  10A40*M(4,3)*B4; 

A2(ab)  =  10A80*M23(1,2); 

B2(ab)  =  10A80*M23(2,2)*B4; 

%  if  ab==l 
%  A1(1)=0; 

%  else 

Al(ab)  =10A20*M123(1,2); 

%  end 
B1  =  0; 

if  F  >  0 

X  =  -300e-10:deltax:0; 

Y  =  0:deltax:L(l,2); 

Z  =  L(l,2):deltax:L(l,3); 

T  =  L(  1 ,3):deltax:(L(  1 ,3)+300e- 1 0); 
else 

X  =  L(  1 ,3):deltax:(L(  1 ,3)+300e- 1 0); 

Y  =  L(l,2):deltax:L(l,3); 

Z  =  0:deltax:L(l,2); 

T  =  -300e-10:deltax:0; 
end 

Pl_bound(ab,:)  =  -C(1).*(F.*X  +  eta(l)); 

P2_bound(ab,:)  =  -C(2).*(F.*Y  +  eta(2)); 

P3_bound(ab,:)  =  -C(3).*(F.*Z  +  eta(3)); 

P4_bound(ab,:)  =  -C(4).*(F.*T  +  eta(4)); 

psil_bound(ab,:)  =  (Al(:,ab).*asymptotieAiry(0,Pl_bound(ab,:))  + 
B 1  ,*asymptoticAiry(2,Pl_bound(ab,:))); 


88 


psi2_bound(ab,:)  =  (A2(:,ab).  ^asymptotic  A  iry(0,P2_bound(ab,:))  + 

B2( :  ,ab) .  *  asymptotic  Airy(2  ,P2_bound(ab, :))); 

psi3_bound(ab,:)  =  (A3(:,ab).*asymptoticAiry(0,P3_bound(ab,:))  + 

B3(:,ab).*asymptoticAiry(2,P3_bound(ab,:))); 

psi4_bound(ab,:)  =  (A4.*asymptoticAiry(0,P4_bound(ab,:))  + 
B4.*asymptoticAiry(2,P4_bound(ab,:))); 

psil_bound(ab,:)  =  10A0*psil_bound(ab,:); 
psi2_bound(ab,:)  =  10A-100*psi2_bound(ab,:); 
psi3_bound(ab,:)  =  10A-100*psi3_bound(ab,:); 
psi4_bound(ab,:)  =  10A-100*psi4_bound(ab,:); 

integrandl_bound(ab,:)  =  (psil_bound(ab,:).*conj(psil_bound(ab,:)))*deltax; 
integrand2_bound(ab,:)  =  (psi2_bound(ab,:).*conj(psi2_bound(ab,:)))*deltax; 
integrand3_bound(ab,:)  =  (psi3_bound(ab,:).*conj(psi3_bound(ab,:)))*deltax; 
integrand4_bound(ab,:)  =  (psi4_bound(ab,:).*conj‘(psi4_bound(ab,:)))*deltax; 
integral  l_bound(ab,:)  =  sum(integrandl_bound(ab,:)); 
integral2_bound(ab,:)  =  sum(integrand2_bound(ab,:)); 
integral3_bound(ab,:)  =  sum(integrand3_bound(ab,:)); 
integral4_bound(ab,:)  =  sum(integrand4_bound(ab,:)); 

finali_bound(ab,:)  =  integral  l_bound(ab,:)  +  integral2_bound(ab,:)  +  inte- 
gral3_bound(ab,:)  +  integral4_bound(ab,:); 

b3_bound(ab,:)  =  l./(sqrt(finali_bound(ab,:))); 

Psil_bound(ab,:)  =  b3_bound(ab,:).*psil_bound(ab,:); 

Psi2_bound(ab,:)  =  b3_bound(ab,:).!|spsi2_bound(ab,:); 

Psi3_bound(ab,:)  =  b3_bound(ab,:).!|spsi3_bound(ab,:); 

Psi4_bound(ab,:)  =  b3_bound(ab,:).*psi4_bound(ab,:); 

Integral  l_bound(ab,:)=sum((Psi  l_bound(ab,:)Aconj(Psi  l_bound(ab,:)))*deltax); 
Integral2_bound(ab,:)=sum((Psi2_bound(ab,:).*conj(Psi2_bound(ab,:)))*deltax); 
Integral3_bound(ab,:)=sum((Psi3_bound(ab,:).*conj(Psi3_bound(ab,:)))*deltax); 
Integral4_bound(ab,:)=sum((Psi4_bound(ab,:).*conj(Psi4_bound(ab,:)))*deltax); 
Ttal_integral_bound(ab,:)=Integrall_bound(ab,:)+Integral2_bound(ab,:)+Integral3_boun 
d(ab,:)+Integral4_bound(ab,:) 


%**Simulation  Of  Performance  Of  Infrared  Photodetectors'1”1”1”1”1”1”1' 
%Epsi_contl 

%LT  Psarakis  Eftychios  Hellenic  Navy 
deltax  =  1 0A- 11; 
if  F>0 

X  =  -300e-10:deltax:0; 

Y  =  0:deltax:L(l,2); 

Z  =  L(l,2):deltax:L(l,3); 

T  =  L(  1 ,3):deltax:(L(  1 ,3)+300e- 1 0); 

A4(ac)  =  l/sqrt(2'1'pi); 

B4(ac)  =  -(M 1 23_c(2, 1  )/M  1 23_c(2,2))* A4(ac)  ; 

A3(ac)  =  (M_c(l,3)*A4(ac)  +  B4(ac)*M_c(2,3)); 
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B3(ac)  =  (M_c(3,3)*A4(ac)  +  B4(ac)*M_c(4,3)); 

A2(ac)  =  (M23_c(  1 , 1  )*A4(ac)+B4(ac)*M23_c(  1 ,2)); 

B2(ac)  =  (M23_c(2, 1  )*  A4(ac)+B4(ac)*M23_c(2,2)); 

Al(ac)  =  (M 1 23_c(  1 , 1  )*  A4(ac)  +  B4(ac)*M123_c(l,2)); 

Bl(ac)  =  0; 
else 

X  =  L(  1 ,3  ):deltax:(L(  1 ,3)+300e- 1 0); 

Y  =  L(l,2):deltax:L(l,3); 

Z  =  0:deltax:L(l,2); 

T  =  -300e-10:deltax:0; 

A4(ac)  =  l/sqrt(2*pi); 

B4(ac)  =  -(M123_c(2,l)/M123_c(2,2))*A4(ac)  ; 

A3(ac)  =  (M_c(l,3)*A4(ac)  +  B4(ac)*M_c(2,3)); 

B3(ac)  =  (M_c(3,3)*A4(ac)  +  B4(ac)*M_c(4,3)); 

A2(ac)  =  (M23_c(  1 , 1  )*A4(ac)+B4(ac)*M23_c(  1 ,2)); 

B2(ac)  =  (M23_c(2, 1  )*  A4(ac)+B4(ac)*M23_c(2,2)); 

Al(ac)  =  (M  123_c(  1 , 1  )*  A4(ac)  +  B4(ac)*M123_c(l,2)); 

Bl(ac)  =  0; 
end 

Pl_cont(ac,:)  =  -C(1).*(F.*X  +  eta(l)); 

P2_cont(ac,:)  =  -C(2).*(F.*Y  +  eta(2)); 

P3_cont(ac,:)  =  -C(3).*(F.*Z  +  eta(3)); 

P4_cont(ac,:)  =  -C(4).*(F.*T  +  eta(4)); 

psil_cont(ac,:)  =  Al(ac)*asymptoticAiry(0,Pl_cont(ac,:))  + 

B 1  (ac)  *  asymptotic  Airy  (2  ,P  1  _cont(  ac , : ; 

psi2_cont(ac,:)  =  A2(ac)*asymptoticAiry(0,P2_cont(ac,:))  + 

B2(ac)*asymptoticAiry(2,P2_cont(ac,:)); 

psi3_cont(ac,:)  =  A3(ac)*asymptoticAiry(0,P3_cont(ac,:))  + 

B3(ac)*asymptoticAiry(2,P3_cont(ac,:)); 

psi4_cont(ac,:)  =  A4(ac)*asymptoticAiry(0,P4_cont(ac,:))  + 

B4(ac)*asymptoticAiry(2,P4_cont(ac,:)); 

integrand  l_cont(ac,:)  =  (psil_cont(ac,:).!|sconj(psil_cont(ac,:)))*deltax; 
integrand2_cont(ac,:)  =  (psi2_cont(ac,:).*conj(psi2_cont(ac,:)))*deltax; 
integrand3_cont(ac,:)  =  (psi3_cont(ac,:).*conj(psi3_cont(ac,:)))*deltax; 
integrand4_cont(ac,:)  =  (psi4_cont(ac,:).*conj(psi4_cont(ac,:)))*deltax; 
integrall_cont(ac,:)  =  sum(integrandl_cont(ac,:)); 
integral2_cont(ac,:)  =  sum(integrand2_cont(ac,:)); 
integral3_cont(ac,:)  =  sum(integrand3_cont(ac,:)); 
integral4_cont(ac,:)  =  sum(integrand4_cont(ac,:)); 

finali_cont(ac,:)  =  integral  l_cont(ac,:)  +  integral2_cont(ac,:)  +  inte- 

gral3_cont(ac,:)  +  integral4_cont(ac,:); 

b3_cont(ac,:)  =  l/(sqrt(finali_cont(ac,:))); 

Psil_cont(ac,:)  =  b3_cont(ac,:)*psil_cont(ac,:); 

Psi2_cont(ac,:)  =  b3_cont(ac,:)*psi2_cont(ac,:); 

Psi3_cont(ac,:)  =  b3_cont(ac,:)*psi3_cont(ac,:); 

Psi4_cont(ac,:)  =  b3_cont(ac,:)*psi4_cont(ac,:); 
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Integral  I_cont(ac,:)=sum((Psi  l_cont(ac,:).*conj(Psi  l_cont(ac,:)))*deltax); 
Integral2_cont(ac,:)=sum((Psi2_cont(ac,:).*conj(Psi2_cont(ac,:)))*deltax); 
Integral3_cont(ac,:)=sum((Psi3_cont(ac,:).*conj(Psi3_cont(ac,:)))*deltax); 
Integral4_cont(ac,:)=sum((Psi4_cont(ac,:).*conj(Psi4_cont(ac,:)))*deltax); 

To- 

tal_integral_cont(ac,:)=Integrall_cont(ac,:)+Integral2_cont(ac,:)+Integral3_cont(ac,:)+Int 

egra!4_cont(ac,:); 


%**Simulation  Of  Performance  Of  Infrared  Photodetectors******* 

%Epsi_bound 

%LT  Psarakis  Eftychios  Hellenic  Navy 
deltax  =  10A-11; 

A4_b  =  0; 

B4_b=  1; 

A3_b(a)  =  M(l,3)*A4_b  +  M(2,3)*B4_b; 

B3_b(a)  =  M(3,3)*A4__b  +  M(4,3)*B4_b; 

A2_b(a)  =  M23(l,2)*B4_b; 

B2_b(a)  =  M23(2,l)*A4_b  +  M23(2,2)*B4_b; 

Al_b(a)  =  M123(l,2)*B4_b; 

Bib  =  0; 
if  F  >  0 

X  =  -300e-10:deltax:0; 

Y  =  0:deltax:L(l,2); 

Z  =  L(l,2):deltax:L(l,3); 

T  =  L(  1 ,3):deltax:(L(  1 ,3)+300e- 1 0); 
else 

X  =  L(l,3):deltax:(L(l,3)+300e-10); 

Y  =  L(l,2):deltax:L(l,3); 

Z  =  0:deltax:L(l,2); 

T  =  -300e-10:deltax:0; 
end 

Pl_bound(a,:)  =  -C(1)A(F.*X  +  eta(l)); 

P2_bound(a,:)  =  -C(2)A(F.sl!Y  +  eta(2)); 

P3_bound(a,:)  =  -C(3)A(F.*Z  +  eta(3)); 

P4_bound(a,:)  =  -C(4).*(F.!|ST  +  eta(4)); 

psil_bound(a,:)  =  Al_b(:,a).*airy(0,Pl_bound(a,:))  + 

B 1  _b .  *  airy  (2  ,P  1  _bound(a, : )) ; 

psi2_bound(a,:)  =  A2_b(:,a).*airy(0,P2_bound(a,:))  + 

B2_b( :  ,a) .  *  airy(2,P2_bound(a, :)); 

psi3_bound(a,:)  =  A3_b(:,a).*airy(0,P3_bound(a,:))  + 

B3_b(:,a).*airy(2,P3_bound(a,:)); 
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psi4_bound(a,:)  =  A4_b.*airy(0,P4_bound(a,:))  +  B4_b.*airy(2,P4_bound(a,:)); 


psi_l_bound(a,:)  =  conj(psil_bound(a,:)); 
psi_2_bound(a,:)  =  conj(psi2_bound(a,:)); 
psi_3_bound(a,:)  =  conj(psi3_bound(a,:)); 
psi_4_bound(a,:)  =  conj(psi4_bound(a,:)); 

integrand l_bound(a,:)  =  (psi  1  _bound(a, : ).  *psi_  1  _bound(a, :))  *  deltax; 
integrand2_bound(a,:)  =  (p  s  i  2_bo  un  d(a, : ) .  *  p  s  i_2_bo  un  d(a, : ) )  *  de  1  tax ; 
integrand3_bound(a,:)  =  (psi3_bound(a,:).*psi_3_bound(a,:))*deltax; 
integrand4_bound(a,:)  =  (psi4_bound(a,:).*psi_4_bound(a,:))*deltax; 

integrall_bound(a,:)  =  sum(integrandl_bound(a,:)); 
integral2_bound(a,:)  =  sum(integrand2_bound(a,:)); 
integral3_bound(a,:)  =  sum(integrand3_bound(a,:)); 
integra!4_bound(a,:)  =  sum(integrand4_bound(a,:)); 


finali_bound(a,:)  =  integral  l_bound(  a,:)  +  integral2_bound(a,:)  +  inte- 

gral3_bound(a,:)  +  integral4_bound(a,:); 


b3_bound(a,:)  =  l/(sqrt(finali_bound(a,:))); 


Psil_bound(a,:)  =  b3_bound(a,:)!|spsil_bound(a,:); 
Psi2_bound(a,:)  =  b3_bound(a,:)*psi2_bound(a,:); 
Psi3_bound(a,:)  =  b3_bound(a,:)*psi3_bound(a,:); 
Psi4_bound(a,:)  =  b3_bound(a,:)!|spsi4_bound(a,:); 


%** Simulation  Of  Performance  Of  Infrared  Photodetectors'1”1”1”1”1”1”1' 
%Epsi_cont 

%LT  Psarakis  Eftychios  Hellenic  Navy 
deltax  =  10A-1 1; 

if  F>0 

X  =  -300e-10:deltax:0; 

Y  =  0:deltax:L(l,2); 

Z  =  L(l,2):deltax:L(l,3); 

T  =  L(  1 ,3):deltax:(L(  1 ,3)+300e- 1 0); 

A4(ac)  =  l/sqrt(2*pi); 

B4(ac)  =  -(M 1 23_c(2, 1  )/M  1 23_c(2,2))* A4(ac)  ; 

A3(ac)  =  (M_c(l,3)!|5A4(ac)  +  B4(ac)'1'M_c(2,3)); 

B3(ac)  =  (M_c(3,3)*A4(ac)  +  B4(ac)*M_c(4,3)); 

A2(ac)  =  (M23_c(l,l)*A4(ac)+B4(ac)*M23_c(l,2)); 


92 


B2(ac)  =  (M23_c(2,l)*A4(ac)+B4(ac)*M23_c(2,2)); 
Al(ac)  =  (M 1 23_c(  1 , 1 )  *  A4(ac)  +  B4(ac)*M123_c(l,2)); 
Bl(ac)  =  0; 
else 

X  =  L(l,3):deltax:(L(l,3)+300e-10); 

Y  =  L(l,2):deltax:L(l,3); 

Z  =  0:deltax:L(l,2); 

T  =  -300e-10:deltax:0; 

A4(ac)  =  l/sqrt(2*pi); 

B4(ac)  =  -(M123_c(2,l)/M123_c(2,2))*A4(ac)  ; 

A3(ac)  =  (M_c(l,3)*A4(ac)  +  B4(ac)*M_c(2,3)); 

B3(ac)  =  (M_c(3,3)*A4(ac)  +  B4(ac)*M_c(4,3)); 

A2(ac)  =  (M23_c(  1 ,  l)*A4(ac)+B4(ac)*M23_c(  1 ,2)); 
B2(ac)  =  (M23_c(2,l)*A4(ac)+B4(ac)*M23_c(2,2)); 
Al(ac)  =  (M123_c(l,l)*A4(ac)  +  B4(ac)*M123_c(l,2)); 
Bl(ac)  =  0; 
end 

Pl_cont(ac,:)  =  -C(1).*(F.*X  +  eta(l)); 

P2_cont(ac,:)  =  -C(2).*(F.*Y  +  eta(2)); 

P3_cont(ac,:)  =  -C(3).*(F.*Z  +  eta(3)); 

P4_cont(ac,:)  =  -C(4).*(F.*T  +  eta(4)); 


psil_cont(ac,:)  =  Al(ac)*airy(0,Pl_cont(ac,:))  +  Bl(ac)*airy(2,Pl_cont(ac,:)); 
psi2_cont(ac,:)  =  A2(ac)*airy(0,P2_cont(ac,:))  +  B2(ac)*airy(2,P2_cont(ac,:)); 
psi3_cont(ac,:)  =  A3(ac)*airy(0,P3_cont(ac,:))  +  B3(ac)*airy(2,P3_cont(ac,:)); 
psi4_cont(ac,:)  =  A4(ac)*airy(0,P4_cont(ac,:))  +  B4(ac)*airy(2,P4_cont(ac,:)); 

psi_l_cont(ac,:)  =  conj(psil_cont(ac,:)); 
psi_2_cont(ac,:)  =  conj(psi2_cont(ac,:)); 
psi_3_cont(ac,:)  =  conj(psi3_cont(ac,:)); 
psi_4_cont(ac,:)  =  conj(psi4_cont(ac,:)); 

integrand  l_cont(ac,:)  =  (psil_cont(ac,:).*psi_l_cont(ac,:))*deltax; 
integrand2_cont(ac,:)  =  (psi2_cont(ac,:).*psi_2_cont(ac,:))*deltax; 
integrand3_cont(ac,:)  =  (psi3_cont(ac,:).*psi_3_cont(ac,:))*deltax; 
integrand4_cont(ac,:)  =  (psi4_cont(ac,:).*psi_4_cont(ac,:))*deltax; 


integrall_cont(ac,:)  =  sum(integrandl_cont(ac,:)); 
integral2_cont(ac,:)  =  sum(integrand2_cont(ac,:)); 
integral3_cont(ac,:)  =  sum(integrand3_cont(ac,:)); 
integral4_cont(ac,:)  =  sum(integrand4_cont(ac,:)); 
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finali_cont(ac,:)  =  integrall_cont(ac,:)  +  integral2_cont(ac,:)  +  inte- 
gral3_cont(ac,:)  +  integral4_cont(ac,:); 

b3_cont(ac,:)  =  l/(sqrt(finali_cont(ac,:))); 

Psil_cont(ac,:)  =  b3_cont(ac,:)*psil_cont(ac,:); 

Psi2_cont(ac,:)  =  b3_cont(ac,:)!|!psi2_cont(ac,:); 

Psi3_cont(ac,:)  =  b3_cont(ac,:)*psi3_cont(ac,:); 

Psi4_cont(ac,:)  =  b3_cont(ac,:)*psi4_cont(ac,:); 

psil_cont_der(ac,:)  =  Al(ac)*airy(l,Pl_cont(ac,:)) 

+B 1  (ac)*airy(3,P  l_cont(ac,:)); 

psi2_cont_der(ac,:)  =  A2(ac)*airy(l,P2_cont(ac,:))  + 

B2(ac)*airy(3,P2_cont(ac,:)); 

psi3_cont_der(ac,:)  =  A3(ac)*airy(l,P3_cont(ac,:))  + 

B3(ac)*airy(3,P3_cont(ac,:)); 

psi4_cont_der(ac,:)  =  A4(ac)*airy(l,P4_cont(ac,:))  + 

B4(ac)*airy(3,P4_cont(ac,:)); 

Psil_cont_der(ac,:)  =  b3_cont(ac,:)*psil_cont_der(ac,:); 

Psi2_cont_der(ac,:)  =  b3_cont(ac,:)*psi2_cont_der(ac,:); 

Psi3_cont_der(ac,:)  =  b3_cont(ac,:)*psi3_cont_der(ac,:); 

Psi4_cont_der(ac,:)  =  b3_cont(ac,:)*psi4_cont_der(ac,:); 

Integral  l_cont(ac,:)=sum((Psi  l_cont(ac,:).*conj(Psi  1  _cont(ac,:)))*deltax); 
Integral2_cont(ac,:)=sum((Psi2_cont(ac,:).*conj(Psi2_cont(ac,:)))*deltax); 
Integral3_cont(ac,:)=sum((Psi3_cont(ac,:).*conj(Psi3_cont(ac,:)))*deltax); 
Integral4_cont(ac,:)=sum((Psi4_cont(ac,:).*conj(Psi4_cont(ac,:)))*deltax); 


To- 

tal_integral_cont(ac,:)=Integrall_cont(ac,:)+Integral2_cont(ac,:)+Integral3_cont(ac,:)+Int 

egra!4_cont(ac,:); 


%**Simulation  Of  Performance  Of  Infrared  Photodetectors******* 
%LT  Psarakis  Eftychios  Hellenic  Navy 
function  y=aproxAiry(x) 

%this  function  aproximates  the  airy's  functions  and  their 
%dirivatives 

%Author:  ©  E. Psarakis,  §  31 -Aug-2004 

%Ai(x)— >ans(l,l)  — >  ai 

%Bi(x)— >ans(l,2)  — >  bi 

%Ai'(x)— >ans(  1 ,3)— >  ad 

%Bi'(x)— >ans(  1 ,4)— >  bd 

%ai=  []  ;bi=[]  ;ad=[]  ;bd=  [] ; 

%ai=(4*pi)A(-l/2)*(x).A(-l/4).*exp(-2/3*(x).A(3/2)); 

%bi=(pi)A(-l/2)*(x).A(-l/4).*exp(2/3!|!(x).A(3/2)); 

%ad=-(4*pi)A(-l/2)*(x).A(l/4).*exp(-2/35|s(x).A(3/2)); 
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%bd=(pi)A(-l/2)*(x).A(l/4).*exp(2/3*(x).A(3/2)); 

%[ai;  bi  ;ad  ;bd]; 

%y=eval('[ai;bi;ad;bd]’); 

index=0; 

ai=[]  ;bi=[]  ;ad=[]  ;bd=[] ; 
for  index=l  :length(x); 
if  x(index)>=0 

ai(index)=(4*pi)A(H/2)*(x(index)).A(-l/4).*exp(-2/3*(x(index)).A1.5); 
bi(index)=(pi)A(-l/2)*(x(index))  A(-l/4).*exp(2/3*(abs(x(index)))  A1.5); 
ad(index)=-(4*pi)A(-l/2)*(x(index)).A(l/4).*exp(-2/3*(x(index)).A1.5); 
bd(index)=(pi)A(-l/2)*(x(index))  .A(l/4).*exp(2/3*(x(index)).A1.5); 
else 

ai(index)=(pi)A(-l/2)*(-(x(index))).A(-l/4).*sin(2/3*(-(x(index)))  A1.5+pi/4); 
bi(index)=(pi)A(-l/2)*(-(x(index))).A(-l/4).*cos(2/3*(-(x(index))).A1.5+pi/4); 
ad(index)=-(pi)A(-l/2)*(-(x(index)))  A(l/4).*(cos(2/3*(-(x(index))).A1.5+pi/4)); 
%-l/(4*sqrt(pi).*(abs(x(index)))A(-5/4)).*sin(2/3*(abs(x(index)))/  1 ,5+pi/4); 

bd(index)=(pi)A(4/2)*(-(x(index)))A(l/4).*(sin(2/3*(-(x(index))).A1.5+pi/4)); 

%-l/(4*sqrt(pi).*(abs(x(index))).A(-5/4)).*cos(2/3*(abs(x(index))).A1.5+pi/4); 

end 

end 

[ai;  bi;ad  ;bd]; 
y=eval('[ai;  bi;ad  ;bd]'); 

%  ai=[];bi=[];ad=[];bd=[]; 


%**Simulation  Of  Performance  Of  Infrared  Photodetectors'1”1”1”1”1”1”15 
%LT  Psarakis  Eftychios  Hellenic  Navy 
function  t=asymptoticAiry(n,x) 

%this  function  aproximates  the  airy’s  functions  and  their 
%dirivatives.  it  is  used  like  the  matlab  airy(x) 

%Author:  ©  E. Psarakis,  §  31-Aug-2004 
%  Ai(x)— >ans(  1,1) 

%Bi(x)— >ans(  1 ,2) 

%Ai'(x)— >ans(  1 ,3) 

%Bi’(x)— >ans(  1 ,4) 
if  n==0 

aproxAiry(x); 
t=ans(l,:); 
end 
if  n==l 

aproxAiry(x); 

t=ans(3,:); 

end 

if  n==2 

aproxAiry(x); 
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t=ans(2,:); 

end 

if  n==3 

aproxAiry(x); 

t=ans(4,:); 

end 

%  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  NP  S  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 
%**  Simulation  Of  Performance  Of  Infrared  Photodetectors******* 

%ZETA 

%LT  Psarakis  Eftychios  Hellenic  Navy 
warning  off 
deltaz  =  10A-11; 

if  F  >  0 

X  =  -300e-10:deltaz:0; 

Y  =  0:deltaz:L(l,2); 

Z  =  L(l,2):deltaz:L(l,3); 

T  =  L(l,3):deltaz:(L(l,3)+300e-10); 
dim=max(T)-min(X); 
else 

X  =  L(l,3):deltaz:(L(l,3)+300e-10); 

Y  =  L(l,2):deltaz:L(l,3); 

Z  =  0:deltaz:L(l,2); 

T  =  -300e-10:deltaz:0; 
dim=max(X)-min(T); 
end 

for  kk=l:ac-l 

Zital(kk,:)=conj(Psil_bound(l,:)).*X.*Psil_cont(kk,:).*deltaz; 

Zita2(kk,:)=conj(Psi2_bound(l,:)).*Y.*Psi2_cont(kk,:).*deltaz; 

Zita3  (kk,  :)=conj  (Psi3_bound(  1 , :)) .  *Z .  *Psi3_cont(kk, :) .  *  deltaz; 

Zita4(kk,  :)=conj  (Psi4_bound(  1  ,:)).*  T .  *Psi4_cont(kk, :) .  *  deltaz; 
OSC_integrall(kk,:)  =  sum(Zital(kk,:)); 

OSC_integral2(kk,:)  =  sum(Zita2(kk,:)); 

OSC_integral3(kk,:)  =  sum(Zita3(kk,:)); 

OSC_integral4(kk,:)  =  sum(Zita4(kk,:)); 

O  S  C_finali(kk, :  )=(0  S  Cintegral  1  (kk, : )+  O  SC_integral2(kk, : )+ 

OSC_integral3(kk,:)+OSC_integral4(kk,:)); 

Zeta_final(kk,:)=(OSC_finab(kk,:)).*conj(OSC_finali(kk,:)); 

%  DOS(kk)=dim./(hbar*pi)*sqrt(m_eff(l)./(2*q*(Eout_cont(kk)- 

min(  V  ( 1 ),  V  (2))))); 

DOS(kk)=dim./(hbar*pi)*sqrt(m_eff(l)./(2*q*(Eout_cont(kk)- 
Eout_cont(  1  )+.00 1 ))); 

end 


96 


APPENDIX  B:  EKF  INITIALIZATION  AND  TRACKING 


In  the  present  Appendix  we  present  the  BOT  initialization  and  tracking  algorithm. 
The  core  of  the  algorithm  is  BOT.m.  Once  this  program  is  executed,  the  loop  INITIALI¬ 
ZATION.  m  is  called,  and  the  user  must  define  the  initialization  bearings. 


***************** ]^j'Y']J^LJ2ATI0N  m 


%******Writen  by  LT  Eftychios  Psarakis  Hellenic  Navy******************* 
%***Thesis  Project:  Simulation  of  Performance  of  Infrared  Photodetectors***** 

%  clc 
%  clear 

Slp=[0;0];  %  Sensor  1  at  origin 

Slp=[10000;0];  %  Sensor  2  at  (10000,0) 

%  bl=input('give  mesured  bearing  1  from  sensor  1  in  degrees:  ');%33 
%  b2=input('give  mesured  bearing  1  from  sensor  2  in  degrees:  ');%  173 
bl_l=bl*2*pi/360;  %convert  to  rad 
bl_2=b2*2*pi/360;  %convert  to  rad 
sigmalb=l*pi/180;  %  1  deg  std  dev  in  bearing 

miss=ab  s(randn*  sigma  1  b) ; 


bl_Sl_plus  =bl_l+miss;  %  bearing  1  from  sensor  1  (+) 
bl_Sl_minus=bl_l-miss;  %  bearing  1  from  sensor  1  (-) 


bl_S2_plus  =bl_2+miss;  %  bearing  1  from  sensor  2  (+) 
bl_S2_minus=bl_2-miss;  %  bearing  1  from  sensor  2  (-) 


%  b_l=input('give  mesured  bearing  2  from  sensor  1  in  degrees:  ');%33.46 
%  b_2=input(’give  mesured  bearing  2  from  sensor  2  in  degrees:  ’);%  1 73 . 1 8 
%  dt=input(’give  time  interval  between  the  two  measurements  :  ’); 
dt=.l; 

b2_l=b_l*2*pi/360;  %convert  to  rad 
b2_2=b_2*2*pi/360;  %convert  to  rad 


b2_Sl_plus=b2_l+miss;  %  bearing  2  from  sensor  1  (+) 
b2_Sl_minus=b2_l-miss;  %  bearing  2  from  sensor  1  (-) 


b2_S2_plus=b2_2+miss;  %  bearing  2  from  sensor  2  (+) 
b2_S2_minus=b2_2-miss;  %  bearing  2  from  sensor  2  (-) 


%Triangle  1:  Bearing  1  Sensor  1(+),  Bearing  1  Sensor  2(+)  using  sine  law 
theta3_l=pi-b  IS  l_plus-(pi-b  l_S2_plus); 
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Dl=((sin(theta3_l)/10000)A-l)*sin(pi-bl_S2_plus);  %distance  from  (0,0) 

x  1  =D  1  *  cos(b  1  _S  1  _plus) ; 
y  1  =D  1  *  sin(b  1  _S  1  _plus); 

%Triangle  1:  Bearing  1  Sensor  1(+),  Bearing  1  Sensor  2(-)  using  sine  law 
theta3_2=pi-b  1_S  l_plus-(pi-b  l_S2_minus); 

D2=((sin(theta3_2)/10000)A-l)*sin(pi-bl_S2_minus);  %distance  from  (0,0) 
x2=D2  *cos(b  1  _S  1  _plus) ; 
y2=D2*sin(b  1_S  l_plus); 

%Triangle  1:  Bearing  1  Sensor  l(-),  Bearing  1  Sensor  2(+)  using  sine  law 
theta3_3=pi-b  IS  l_minus-(pi-b  l_S2_plus); 

D3=((sin(theta3_3)/10000)A-l)*sin(pi-bl_S2_plus);  %distance  from  (0,0) 
x3=D3  *cos(b  IS  l_minus); 
y3=D3*sin(bl_S  l_minus); 

%Triangle  1:  Bearing  1  Sensor  l(-),  Bearing  1  Sensor  2(-)  using  sine  law 
theta3_4=pi-b  IS  l_minus-(pi-b  l_S2_minus); 

D4=((sin(theta3_4)/10000)A-l)*sin(pi-bl_S2_minus);  %distance  from  (0,0) 
x4=D4*cos(b  1_S  l_minus); 
y4=D4*sin(bl_S  I  _minus) ; 


O/q  clisp^*  ^  ^  ^  ^  ^  ^  jDctrctllGlo^rciiYi^'  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ * 

posl=[xl;yl]; 

pos2=[x2;y2]; 

pos3=[x3;y3]; 

pos4=[x4;y4]; 

%  disp('*******************************************'); 
POSl=[(xl+x2+x3+x4)/4; 

(y  1  +y2+y3+y4)/4] ;  %  middle  of  the  rectungular  shape 


d=max(abs((x2-x3)/(atan((y2-y3)/(x2-x3)))),abs((x4-xl)/(atan((y4-yl)/(x4- 

xl ))))); 


xx=linspace(  1 , 1 0000, 100); 

yy  1  =tan(b  1  _S  1  _plus)*xx; 
yy2=tan(b  1_S  l_minus)*xx; 


a3=y3/(x3- 10000); 

b3=-10000*a3; 

yy3=a3*xx+b3; 
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a4=y4/(x4- 10000); 

b4=-10000*a4; 

yy4=a4*xx+b4; 

%  figure 

%  plot(xx,yyl,'r') 

%  hold  on 
%  plot(xx,yy2,’b’) 

%  hold  on 
%  plot(xx,yy3,'r') 

%  hold  on 
%  plot(xx,yy4,'b') 

%Triangle  2:  Bearing  1  Sensor  1(+),  Bearing  1  Sensor  2(+)  using  sine  law 
theta3_  1  _2=pi-b2_S  1  _plus-(pi-b2_S2_plus); 

Dl_2=((sin(theta3_l_2)/10000)A-l)*sin(pi-b2_S2_plus);  %distance  from  (0,0) 
xl_2=Dl_2*cos(b2_S  l_plus); 
y  1  _2=D  1  _2  *  sin(b2_S  l_plus); 


%Triangle  2:  Bearing  1  Sensor  1(+),  Bearing  1  Sensor  2(-)  using  sine  law 
theta3_2_2=pi-b2_Sl_plus-(pi-b2_S2_minus); 

D2_2=((sin(theta3_2_2)/10000)A-l)*sin(pi-b2_S2_minus);  %distance  from  (0,0) 
x2_2=D2_2*cos(b2_S  l_plus); 
y2_2=D2_2*sin(b2_S  l_plus); 

%Triangle  2:  Bearing  1  Sensor  l(-),  Bearing  1  Sensor  2(+)  using  sine  law 
theta3_3_2=pi-b2_Sl_minus-(pi-b2_S2_plus); 

D3_2=((sin(theta3_3_2)/10000)A-l)*sin(pi-b2_S2_plus);  %distance  from  (0,0) 
x3_2=D3_2*cos(b2_S  l_minus); 
y3_2=D3_2  *  sin(b2_S  l_minus); 

%Triangle  2:  Bearing  1  Sensor  l(-),  Bearing  1  Sensor  2(-)  using  sine  law 
theta3_4_2=pi-b2_Sl_minus-(pi-b2_S2_minus); 

D4_2=((sin(theta3_4_2)/10000)A-l)*sin(pi-b2_S2_minus);  %distance  from  (0,0) 
x4_2=D4_2*cos(b2_S  l_minus); 
y  4_2=D4_2  *  sin(b2_S  l_minus); 


o/'q  dispC  *  ********  second  pdrdllclo^rdin^ 

posl_2=[xl_2;yl_2]; 

pos2_2=[x2_2;y2_2]; 

pos3_2=[x3_2;y3_2]; 

pos4_2=[x4_2;y4_2]; 

%  disp('*******************************************'); 


POS2=[(xl_2+x2_2+x3_2+x4_2)/4; 
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(yl_2+y2_2+y3_2+y4_2)/4];  %  middle  of  the  rectungular  shape 
dl=max(abs((x2_2-x3_2)/(atan((y2_2-y3_2)/(x2_2-x3_2)))),abs((x4_2- 
x  1  _2)/(atan((y4_2-y  1  _2 )/ (x4_2-x  1  _2))))); 


sigma  1  r=max(d,d  1 ); 

xx=linspace(  1 , 1 0000, 100); 

yy  1  _2=tan(b2_S  1  _plus)*xx; 
yy2_2=tan(b2_S  l_minus)*xx; 


a3_2=y3_2/(x3_2- 10000); 

b3_2=-10000*a3_2; 

yy3_2=a3_2*xx+b3_2; 

a4_2=y4_2/(x4_2 - 1 0000); 
b4_2=-10000*a4_2; 
yy4_2=a4_2 *xx+b4_2 ; 

%  figure 

%  plot(xx,yyl_2,'r') 

%  hold  on 
%  plot(xx,yy2_2,’b’) 

%  hold  on 
%  plot(xx,yy3_2,'r') 

%  hold  on 
%  plot(xx,yy4_2,’b’) 

xi=[POS2(  1 , 1  );POS2(2,l);POS  1  ( 1 ,  l);POS  1(2,1)]; 

% 

%  figure 

%  plot(xx,yyl,'r') 

%  hold  on 
%  plot(xx,yy2,’r’) 

%  hold  on 
%  plot(xx,yy3,'r') 

%  hold  on 
%  plot(xx,yy4,'r') 

%  hold  on 
%  plot(xx,yyl_2,'b') 

%  hold  on 
%  plot(xx,yy2_2,'b') 

%  hold  on 
%  plot(xx,yy3_2,’b’) 

%  hold  on 
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%  plot(xx,yy4_2,’b') 

% 

%  text(POS  1  ( 1 , 1  ),POS  1  (2, 1  ),’POS  1  ’) 
%  text(P0S2(  1 , 1  ),POS2(2, 1  ),’P0S2’) 


o^^H^********** * Writen  by  LT  Eftychios  Psarakis  Hellenic  Navy********** 
%***TheSiS  Project:  Simulation  of  Performance  of  Infrared  Photodetectors*** 

clear, clc 

disp('for  no  mesurement  noise  select  O') 
disp('for  mesurement  noise  select  1') 

mnoiseflag=input('SELECT  NOISE  OR  NO  NOISE  CASE:’); 

disp('for  single  sensor  tracking  select  O’) 
disp('for  two  sensor  tracking  select  1') 

twosensorflag=input('SELECT  ONE  OR  TWO  SENSORS  CASE:’); 

s2p  =  [10000;0];  %  Sensor  2  position  (Sensor  1  is  at  (0/0) 

%rand(’seed',0);  , 
delta  =  0.1;  %Time 
nsamples  =  1000; 
initialization; 
x  =  xi; 
q=  1000; 

Q=  [deltaA3/3,  deltaA2/2;  deltaA2/2,  delta]; 


Qi  =  [Q,  zeros(2);  zeros(2),Q]; 

Q  =  q*Qi; 

sb  =  1  *pi/l  80; 
D=[l, 

o, 

o, 

0; 

1/dt, 

o, 

-1/dt, 

0; 

o, 

1, 

o, 

0; 

o, 

x_sp=D*xi; 

1/dt, 

0, 

-1/dt]; 

speed=sqrt((x_sp(2,l))A2+(x_sp(4,l))A2)*2000/3600 
wl=input(’SELECT  THE  DESIRED  ACCELERATION  IN  G’) 
w  =  wl*  10/speed; 

R  =  sb; 
xout  =  [] ; 
tout  =  []; 
terr  =  []; 
chi2err  =  []; 

%  xhat  =  x  +  [100*randn;5*randn;100*randn;5*randn]; 
xhat  =  x  +  0.6*[100;5;100;5]; 
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phat=[dlA2  ,0  ,0  ,0; 

0  ,dlA2 ,0  ,0; 

0  ,0  ,dA2  ,0; 

0  ,0  ,0  ,dA2] 


F=[0,  1,  0, 

0,  0,  0 

0,  0,  0, 

0,  w,  0, 

Fturn  =  expm(F*  delta) 

F  =[1,  delta,  0, 

0,  1,  0 

0,  0,  1, 

0,  0,  0, 


0; 

-w; 

i; 

0]; 


0; 

0; 

delta; 

i]; 


Hp  =  [1  0  0  0;  0  0  1  0]; 
for  ii  =l:nsamples 
%  Take  a  measurement 

b  =  atan2(x(3),x(l)); 
if  mnoiseflag==  1 
z  =  b  +  sb*randn; 
else 
z  =  b; 
end 

xout  =  [xout,  [x(l)  ;x(3)  ]  ]; 

%  EKF  Update 
r2  =xhat(l)A2  +  xhat(3)A2; 
b  =  atan2(xhat(3),xhat(l)); 
zhat  =  b; 

H=  [-xhat(3)/r2,  0,  xhat(l)/r2,0]; 
Pzi=inv(H*phat*H'  +  R); 
zt  =  z  -  zhat; 
chi2  =zt' *Pzi*zt; 

K  =  phat*H'*Pzi; 

K2  =  eye(4)-K*H; 
xhat  =xhat  +  K*zt; 
phat  =  K2*phat*K2'  +  K*R*K’; 


if  twosensorflag  ==1  %  process  measuransnt  from  second  sensor 

%  Take  a  Measuronent 
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rp  =Hp*x  -  s2p; 

b  =  atan2(rp(2),rp(l)); 

if  mnoiseflag  ==  1 

z  =  b  +  sb*randn; 

else 

z  =  b; 

end 

%  Measuranent  update 
rp  =  Hp*xhat  -  s2p; 
r2  =  rp'*rp; 
b  =  atan2(rp(2),rp(l)); 
zhat  =  b; 

H=  [-rp(2)/r2,  0,rp(l)/r2,0]  ; 

Pzi  =  inv(H*phat*H’  +  R)  ; 
zt  =  z-  zhat; 
chi2  =  zf  *Pzi*zt; 

K  =  phat*H'*Pzi; 

K2  =  eye(4)  -  K*H; 
xhat  =  xhat  +  K*zt; 
phat  =  K2*phat*K2'  +  K*R*K'; 


end 

e  =  [xhat(l);xhat(3)]; 
tout  =[tout,e] ; 
e  =  e  -  [x(l);x(3)]; 
terr=  [terr,sqrt(e'*e)]; 
chi2err  =  [chi2err,chi2]; 

%  Target  Motion 

if  ((ii  >  250)  &  (ii  <=  350)) 

x  =  Fturn*x; 

else 

x  =  F*x; 
end 

%  Track  Prediction 
xhat=F*xhat; 
phat=F*phat*F'  +  Q; 
end  %  for  ii 

if  mnoiseflag==0  &&  twosensorflag==0 

tlabel  =  ['BO  Tracking  single  Sensor, no  Noise,qA2  =  1000,  accel 

’,num2str(wl),'g']; 

elseif  mnoiseflag==l  &&  twosensorflag==0 

tlabel  =  [’B0  Tracking  single  Sensor, Noise  presence, qA2  =  1000,  accel 
’,num2str(wl),’g']; 

elseif  mnoiseflag==0  &&  twosensorflag==l 

tlabel  =  ['B0  Tracking  two  Sensors,  no  Noise,  qA2  =  1000,  accel 

’,num2str(wl),'g']; 
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elseif  mnoiseflag==l  &&  twosensorflag==l 

tlabel  =  ['BO  Tracking  two  Sensors, Noise  presence,qA2  =  1000,  accel 
’,num2str(wl),'g']; 
end 
figure 

subplot(2,l,l) 

plot(xout(l,:),xout(2,:),  ’:’,tout(l,  :),tout(2,:),’r’); 

%title(  ’Target  and  Track’)  ; 
title  ([tlabel]) ; 
xlabel(’meters’) 
ylabel(’meters’) 

legend(’true  position’, 'Track  position’) 
grid 

subplot(2,l,2) 

plot(xout(l,:),xout(2,:),  ",tout(l,  :),tout(2,:),’r'); 

%title(  ’Target  and  Track’)  ; 
title  ([tlabel]) ; 
xlabel(’meters’) 
ylabel(’meters’) 

legend(’true  position’, ’Track  position’) 
axis([.6*  1 0A4,3 .3  *  1 0A4,3  *  1 0A4, 13*1 0A4]) 
grid 
figure 

tt  =  [0:  (max  (size(terr))-l)]  *delta; 

subplot(2,l,l) 

plot  (tt,terr,  ’r’); 

title  ([  ’Track  Position  Errors:  ’, tlabel]) ; 
xlabel(’Sample  number’) 
ylabel(’meters’) 
grid 

subplot(2,l,2) 

tt  =  [0:  (max(size  (chi2err))-l)]  *delta; 
plot  (tt,chi2err,  ’r’)  ; 

title  (’Chi-Squared  Values  for  Measurement  Association’); 

xlabel(’Sample  number’) 

ylabel('rad') 

grid 
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APPENDIX  C:  ASYMPTOTIC  EXPRESSIONS 


Appendix  C  presents  the  Airy  function  considerations  for  the  low  bias  simulation 
of  the  quantum  well. 

A.  AIRY  FUNCTION  ANALYSIS 

1.  General 


The  differential  equation — --zco  =  0  ,  where  Q  =  2/3  z3/2  and  co  =  e2mn  ,  can  be 

dz 

reduced  to  the  differential  equation  satisfied  by  Bessel  functions  of  order  1/3  [33,  sec. 
6.4].  There  are  two  linearly  independent  solutions  which  are  called  Airy  functions  of  the 
first  and  the  second  kind. 


Ai(z)  =  „(£)-!,„(£)] 

(C.l) 

a'(z)=(/%)  U-m(()+ Wf)] 

(C.2) 

(C.3) 

Airy  Unctions  of  the  first  kind 

(C.4) 

-0.6 


-0.8 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

-20  -15  -10  -5  0  5  10  15  20 

x 

Figure  C.  1 :  Airy  functions  of  the  first  kind,  Ai(-z)  and  Ai(z ) 
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Airy  functions  of  the  second  kind 


Figure  C.2:  Airy  functions  of  the  second  kind,  Bi(-z)  and  Bi(z) 


2.  Asymptotic  Expressions  of  Airy  Functions 


There  are  several  problems  associated  with  the  direct  approach  to  the  solution  of 
the  Schrodinger  equation  using  Airy’s  functions.  The  main  drawback  is  the  asymptotic 
behavior  of  Airy’s  functions  and  of  their  derivatives,  when  the  used  argument,  is  either 
very  large  or  very  small.  In  those  cases  personal  computers  give  numerical  overflows 
which  significantly  reduce  the  result  quality  and  computational  efficiency.  These  kinds  of 
problems  appear  especially  when  low  biases  are  applied  in  symmetric  or  asymmetric 
quantum  wells  structures. 

The  asymptotic  forms  of  Airy’s  functions  used  for  our  simulations  were  taken 
from  the  handbook  of  mathematical  functions  [28,  pp.  448-449]. 

-i  _L  -2  J 

Ai(z)  =  (4  ft)  2z  4e  3  ,z»0  (C.5) 


-i  i  -2zi 

Ai(z)  =  -( 4  ft)  2z4e  3  ,z»0 


Ai(z )  =  (ft)  2  (-z)  4  sin 


( 0  3  A 

2  y  ft 

-(-z)2  +  — 

2  a 

V J  V 


,z  «  0 


(C.6) 

(C.7) 
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Ai  (z)  =  -(a)  2 ( — z)4 cos 


f  ~  3  A 

2  .  n 

-(- z )2  +  — 

3  4 

v J  v 


,z  «  0 


-1  -1  2,:- 

Bi(z)  =  (n )  2z  4e3  ,z»0 

a  i 2  J 

5/  (z)  =  (;r)  2z4e3  ,z»0 


_i  _y  A 

Bi(z )  =  (;r)  2  (-z)  4  cos 


o  3  A 

2  ,  -  T  ^ 

-(-z)2  +  — 

3  4 


_1  1 

Bi  (z)  =  (X)  2  (-z)4  sin 


At  3  A 

h-*r+- 

v3  4 


,z  «  0 


,z  «  0 


(C.8) 

(C-9) 

(C.10) 

(C.ll) 

(C.  12) 


The  expressions  above  give  acceptable  results  not  only  when  the  argument  is  very 
large  or  extremely  small  but  also  when  the  argument  is  a  finite  number  (not  around  zero) 
as  shown  in  Figures  C.3  and  C.4.  In  the  simulation  code,  written  in  Matlab,  we  have  used 
a  specific  procedure  so  that  when  the  built-in  Matlab  Airy  functions  give  results,  we  take 
advantage  and  use  them,  but  when  the  built-in  functions  are  unable  to  give  results  we  use 
Equations  (C.5)  through  (C.12). 


Matlab  Airy  Vs  Asymptotic  Expressions  of  the  first  kind 


1  1  1  1  1  1 

1  1  1  1  -  Matlab  Ary 

|  -  Asymptotic  Expressions 

i  i  i 

i  i  i 

1 

1 

1 

al/ALy . 

v  'O' 

_ 1 _ [ _ 1 _ 

1 

_ 1 _ 

_ 

5 1 _ i _ i _ _ _ i _ i _ i _ i _ i _ 

-8-6-4  -2  02  46  8 


x 


Matlab  Ary  Derivative  Vs  Asymptotic  Expressions  of  the  first  kind 


Figure  C.3:  Matlab  Airy  functions  vs.  asymptotic  expressions  of  the  first  kind 
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15 


Matlab  Airy  Vs  Asymptotic  Expressions  of  the  second  kind 


Matlab  Ary  Derivative  Vs  Asymptotic  Expressions  of  the  second  kind 


Figure  C.4:  Matlab  Airy  functions  vs.  asymptotic  expressions  of  the  second  kind 
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APPENDIX  D:  BOT  COMPLETE  CASE  SET 


We  present  the  complete  results  of  the  BOT  simulation  for  0  g,  2  g,  4  g,  6  g  and  8 
g  turns,  using  one  or  two  sensors,  with  and  without  the  presence  of  measurement  noise. 
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j4  BO  Tracking  single  Sensor,no  Noise, q2  =  1000,  accel  =2g 
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